close

Вход

Забыли?

вход по аккаунту

?

Microwave heating of fluid/solid layers: A study of hydrodynamic stability and melting front propagation

код для вставкиСкачать
INFORMATION TO USERS
This manuscript has been reproduced from the microfilm master. UMI
films the text directly from the original or copy submitted. Thus, some
thesis and dissertation copies are in typewriter face, while others may be
from any type o f computer printer.
T he quality o f this reproduction is dependent upon the quality o f the
copy subm itted.
Broken or indistinct print, colored or poor quality
illustrations and photographs, print bleedthrough, substandard margins,
and improper alignment can adversely affect reproduction.
In the unlikely event that the author did not send UMI a complete
manuscript and there are missing pages, these will be noted.
Also, if
unauthorized copyright material had to be removed, a note will indicate
the deletion.
Oversize materials (e.g., maps, drawings, charts) are reproduced by
sectioning the original, beginning at the upper left-hand comer and
continuing from left to right in equal sections with small overlaps. Each
original is also photographed in one exposure and is included in reduced
form at the back o f the book.
Photographs included in the original manuscript have been reproduced
xerographically in this copy. Higher quality 6” x 9” black and white
photographic prints are available for any photographs or illustrations
appearing in this copy for an additional charge. Contact UMI directly to
order.
UMI
A Bell & Howell Information Company
300 North Zeeb Road, Ann Arbor MI 48106-1346 USA
313/761-4700 800/521-0600
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
NOTE TO USERS
The original manuscript received by UMI contains pages with
slanted print. Pages were microfilmed as received.
This reproduction is the best copy available
UMI
V
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
M IC R O W A V E H E A T IN G O F F L U ID /S O L ID L A Y E R S: A S T U D Y O F
H Y D R O D Y N A M IC S T A B IL IT Y A N D M E L T IN G F R O N T
P R O P A G A T IO N
by
John G ilchrist
A D issertation
Subm itted to th e Faculty o f
N ew Jersey In stitu te o f T echnology
in P artial Fulfillment o f the R equirem ents for th e D egree o f
D octor o f Philosophy
D epartm ent o f M athem atics
A ugust 1998
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
UMI Number: 9906988
Copyright 1998 by
Gilchrist, John Joseph
All rights reserved.
UMI Microform 9906988
Copyright 1998, by UMI Company. All rights reserved.
This microform edition is protected against unauthorized
copying under Title 17, United States Code.
UMI
300 North Zeeb Road
Ann Arbor, M l 48103
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Copyright ©
1998 by John Gilchrist
ALL RIGHTS RESERVED
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
APPRO VAL PAGE
M IC R O W A V E H E A T IN G O F F L U I D /S O L I D L A Y E R S : A S T U D Y O F
H Y D R O D Y N A M I C S T A B IL IT Y A N D M E L T IN G F R O N T
P R O P A G A T IO N
John G ilchrist
-------— ___________________ Q y*U
fir
G regory A . K riegsm ann, Ph.D ., D issertation Advisor
D istinguished Professor, D epartm ent of M athem atics,
New Jersey In stitu te of Technology, Newark N J
2'2i
D ate
*./ I t / 17 1
—DelHeEriuTPapageorgiou, Ph.D ., D issertation Advisor
A ssociate Professor, D epartm ent of M athem atics.
New Jersey In stitu te of Technology, Newark N J
%
D ate
Jon#£fi£in H. C. Luke, Ph.D ., C om m ittee M ember
A ssociate Professor, D epartm ent of M athem atics,
New Jersey In stitu te of Technology, Newark N J
D ate
Yuriko Renardy, Ph.D ., Com m ittee M em ber
Professor, D epartm ent of M athem atics,
V irginia Polytechnic In stitu te and S tate University, Blacksburg VA
D ate
Z?A}4L ^6
B uirt
rt S. Tilley, Ph.D ., Commi£te<
CommiJ^ee M em ber
A ssistant Professor, D epartm ent of M athem atics.
New Jersey In stitu te of Technology, Newark NJ
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
D ate
A B ST R A C T
M IC R O W A V E H E A T IN G O F F L U ID /S O L ID L A Y E R S: A S T U D Y OF
H Y D R O D Y N A M IC S T A B IL IT Y A N D M E L T IN G F R O N T
P R O P A G A T IO N
by
John Gilchrist
In this work we study th e effects of externally induced heating on th e dynamics
of fluid layers, and m aterials composed of two phases separated by a therm ally driven
moving front. One novel aspect of our study, is in the nature of the external source
which is provided by th e action of microwaves acting on dielectric m aterials. The
main challenge is to model and solve systems of differential equations which couple
fluid dynam ical motions (the Navier-Stokes equations for non-isotherm al flows) and
electromagnetic wave propagation (governed by Maxwell’s equations).
W hen an electromagnetic wave impinges on a m aterial, energy is generated
w ithin the m aterial due to dipolar and ohmic heating. The electrical and therm al
properties of the m aterial dictate the dynamics of the heating process, as well
as steady-state tem perature profiles.
Such forms of heating have received little
attention in studies of hydrodynam ic instabilities of non-isotherm al flows, such as the
classical Benard problem, for instance. The novel feature, which allows possibilities
for fluid m anagement and control, is the non-local coupling between th e electro­
m agnetic field and the tem perature distribution within th e fluid. In th e first part
of the thesis, we consider hydrodynamic instabilities of such systems w ith particular
emphasis on conditions for onset of convection. This is achieved by solving th e linear
stability equations in order to identify param eter values which produce instability.
The analysis and subsequent numerical solutions are carried o ut both for m aterials
w ith constant dielectric attributes (in such cases the electric field equations decouple
and they can be solved in closed form), and m aterials w ith tem perature dependent
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
conductivities, dielectric perm ittivities and dielectric loss factors. In th e latter case
we incorporate known d ata for w ater or ethanol into our numerical solutions. O ur
solutions provide a complete picture of onset conditions as a function of in p u t power
levels and microwave frequency (or equivalently fluid layer thickness). In addition,
in the case of water, the flow is found to be more stable for constant attrib u tes as
compared with tem perature dependent attributes; th a t is, a higher power is required
to set the fluid layer into convective motions in the la tte r case. We have also estab­
lished th a t onset is obtained a t power levels well below those needed to cause therm al
runaway and consequently boiling of the w ater layer, for instance. O ur results also
identify different param eter ranges which can produce convection cells of different
sizes w ith the same power input. Such results are directly related to th e microwave
radiation which provides the heating and in particular the distribution of the electric
field w ithin the fluid layer. Several interesting experiments are suggested by the
theoretical predictions.
The second problem studied is concerned with th e use of microwave radiation
in the processing of materials which contain two phases separated by a moving
front which forms and propagates due to a jum p in tem perature flux across the
interface separating the phases. The problem is an extension of th e classical Stefan
problem with the propagation caused by tem perature gradients induced by the
electromagnetic radiation. We have modeled and solved th e problem of two phases
separated by a planar interface and in th e absence of fluid m otion if m elting is
involved.
The boundary conditions are those of convective cooling a t the top
surface and either a heat sink (to m aintain a frozen state for ice, for example) or
a perfectly electromagnetically reflective bounding surface at th e bottom . Known
d a ta modeling a water-ice system have been used, b u t th e m ethods are the same
for other materials.
We have addressed the cases of constant and tem perature
dependent m aterial dielectric attributes.
Steady state configurations have been
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
obtained by solving a coupled system of nonlinear differential equations leading to
an eigenvalue problem for th e interfacial position. In addition, a tim e-dependent
code was developed in order to study transient motions towards steady-state, startin g
from initial configurations of a thin water layer on the ice, for example. O ur results
indicate th a t for a given power level there can be two stable steady-state positions
for the m elting front as well as an unstable one. Existence of m ultiple states is a
consequence of electromagnetic wave resonances within the m aterial and their global
effects on the therm al distribution. Such behavior leads to a theoretical framework
in efforts to control the position of phase separation interfaces in th e processing of
m aterials.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
B IO G R A PH IC A L SK E T C H
A uthor:
John Gilchrist
Degree:
D octor of Philosophy
D ate:
A ugust 1998
D a te o f B ir t h :
July 28, 1963
P la c e o f B ir th :
Kearny, New Jersey
U n d e r g r a d u a te a n d G r a d u a te E d u c a tio n :
• D octor of Philosophy in Applied M athematics,
New Jersey Institute of Technology, Newark, N J, 1998
• M aster of Science in Applied M athematics,
New Jersey In stitu te of Technology, Newark, NJ, 1991
• B atchelor of Science in Computer Science,
R utgers University (Cook College), New Brunswick, N J, 198G
M ajor: Applied Mathematics
P ublications:
John J. Gilchrist, Gregory A. Kriegsmann, Demetrius Papageorgiou, “Stability of
A Microwave H eated Fluid Layer” , IM A Journal o f Applied Mathematics, January
1998.
John J. Gilchrist, C. Musante, S. Nanda, J. Rodricks, D. Trubatch, “Modeling
of T hin Filam ents of Polymeric Liquid Crystals”, Claremont Colleges Mathematics
Modeling Workshop Proceedings, Claremont, CA, June, 1994.
iv
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
To my parents
v
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ACKNOW LEDGM ENT
I, in all probability, never would have pursued a doctoral degree w ithout the
inspiration and encouragement of two very special people, Professor D em etrius
Papageorgiou and Professor Gregory A. Kriegsmann, who saw in me talents and
abilities th a t I didn’t see in myself. I would like to th an k Professor Jonathan Luke,
Professor Yuriko Renardy and Professor B urt S. Tilley who where always there to
offer me guidence in my research. T heir words were always words of encouragement
and these words provided me with th e strength to endure the tough times. I would
like to thank my family who offered their constant support and encouragement.
T hanks to my friends who helped me to m aintain a proper balance between work
and play. They helped me to m aintain my sanity during these years of intense work.
vi
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
TABLE OF C O N T E N T S
C hapter
Page
IN T R O D U C T IO N .........................................................................................................
1
1 .1
Hydrodynamic S tab ility ...................................................................................
1
1.2
Applications of Microwaves in the Melting of M aterials.........................
4
2 HYDRODYNAMIC STABILITY O F A FLUID LAYER: CONSTANT
D IELECTRIC PERM ITTIV ITY C A S E .........................................................
6
1
2.1
The M o d e l..........................................................................................................
2.2
Basic States
.....................................................................................................
10
A Limiting Case k —>0 (Low F requency)........................................
11
2.3
Linear Stability ................................................................................................
13
2.4
Numerical S o lu tio n s ........................................................................................
14
2.2.1
3
6
2.4.1
A Limiting Case (Low Wave N um ber P erturbative Modes) . . .
16
2.4.2
A Limiting Case (Stability of Modesfor Large Critical Power
L e v e ls ).................................................................................................
17
2.5
R e s u l t s ...............................................................................................................
19
2.6
D isc u ssio n ..........................................................................................................
26
HYDRODYNAMIC STABILITY O F A FLUID LAYER: TEM PER A TU R E
D EPEN D EN T COMPLEX P E R M IT T IV IT Y ...............................................
28
3.1
The M o d e l..........................................................................................................
28
3.2
Basic S t a t e s ........................................................................................................
31
3.2.1
Therm al R unaw ay..................................................................................
34
3.3
Linear Stability T h e o ry ...................................................................................
39
3.4
Numerical S o lu tio n s ........................................................................................
42
3.4.1
A Finite Difference A p p ro ach .............................................................
42
3.5
R e s u l t s ................................................................................................................
45
3.6
D isc u ssio n ..........................................................................................................
50
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Chapter
Page
4 A MODEL F O R MELTING OF SOLIDS USING M IC R O W A V E S...............
52
4.1
A Physical Overview (The Stefan C o n d itio n )...........................................
52
4.2
The M athem atical Model .............................................................................
54
4.3
The Case of C onstant Dielectric P erm ittiv ity ...........................................
57
4.4
4.3.1
Steady-State Solutions .......................................................................
57
4.3.2
A sym ptotic Limit of Large Stefan Number ..................................
62
The Case of Tem perature Dependent Dielectric A ttr ib u te s .................
66
4.4.1
4.5
Steady-State Solutions.........................................................................
66
A Fixed Front Method For Tracking the Moving B o u n d ary .................
68
4.5.1
Numerical Im plem entation.................................................................
69
4.5.2
The Case of Constant Dielectric Properties (General Stefan
N u m b e r)..............................................................................................
70
The Case of Tem perature Dependent Dielectric Properties
(General Stefan N u m b e r)...............................................................
72
D isc u ssio n ..........................................................................................................
77
4.5.3
4.6
APPEND IX A
THE M ETHOD OF MULTIPLE SC A LES...................................
78
APPEND IX B
DERIVATION OF EQUATIONS GOVERNING “COM PLEX
ELECTRIC FIELD” ....................................................................
81
WKB A N A L Y S IS..............................................................................
83
REFERENCES ................................................................................................................
85
APPEND IX C
viii
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
L IST OF F IG U R E S
Figure
P age
2.1
The geom etry of the hydrodynamic problem .....................................................
7
2.2
Typical undisturbed tem perature profiles; the physical param eters are
P = 1; er = 2.54 ; f = .25...................................................................................
12
Growth rate curves, rigid-rigid case; P = 1; k = 1; eT = 2.54; P r = 17;
f = .25.....................................................................................................................
15
Eigenfunctions W (z) (top) and 0(.z) (bottom ); P = 1; k = 1; er = 2.54;
P r = 17; f = .25...................................................................................................
17
2.5
A sym ptotic eigenfunctions a t large i?x; P =
19
2.6
N eutral stability curves in the i?x_a plane; k = 1; er = 2.54; ; f = .25. . .
21
2.7
N eutral curves rigid-rigid case: p = 1; er = 2.54; £ = .25...............................
21
2.8
As in Figure 2.7 but rigid-free boundary conditions........................................
22
2.9
P erturbation flow field at onset; base tem perature profile x 10 super­
imposed on flow field. /? = 1; k = .705; er = 2.54; £ = .25; critical
wave num ber ac =2.45........................................................................................
24
2.10 Perturbation flow field at onset; base tem perature profile x 10 super­
imposed on flow field, ft = 1; k = 1.355; er = 2.54; £ = .25; critical
wave number ac =4.401......................................................................................
24
2.11 Perturbation flow field at onset; base tem perature profile x 10 super­
imposed on flow field. /? = 1; k = .71; er = 2.54; £ = .25; critical wave
num ber ac =1.99...................................................................................................
25
2.12 P erturbation flow field at onset; base tem perature profile x 10 super­
imposed on flow field. P = 1; k = 1.36; er = 2.54; £ = .25; critical
wave num ber ac =3.22........................................................................................
25
3.1
Dielectric perm ittivity of water vs. te m p e ra tu re ;...........................................
33
3.2
Dielectric loss factor of water vs. tem perature; ..............................................
33
3.3
U ndisturbed tem perature profiles for water; /3 = 1............................................
34
3.4
Steady-state curves for water (8m vs. x); P = 1................................................
36
3.5
Steady-state curves for w ater (do vs. Xrun); P «
1 ......................................
36
3.6
Growth rate curves, rigid-rigid case; p = 1; k = 1; P r = 7............................
45
2.3
2.4
1;
k =
1;
er = 2.54; f = .25.
ix
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
F ig u r e
Page
3.7
N eutral stability curves in the x~a plane for water; k = 1..............................
47
3.8
N eutral curves in the X c~ k plane for water, rigid-rigid case: /? = 1; k = 1.
48
3.9
N eutral curves for w ater (critical wave num ber vs. k), rigid-rigid case:
/? = 1 ; k = 1...........................................................................................................
49
3.10 N eutral curves for w ater (height of m axim um tem perature vs. k), rigidrigid case: /? = 1; k = 1.......................................................................................
49
3.11 N eutral curves for w ater (maximum tem perature vs. k), rigid-rigid case:
/3 = 1; k = 1...........................................................................................................
50
4.1
The geom etry of the melting problem ..................................................................
53
4.2
The physics governing th e moving front..............................................................
53
4.3
Steady-state curves (constant complex perm ittivity w ith varying Biot
number): Tbw = —.5; k = 1................................................................................
60
Steady-state curves (constant complex perm ittivity w ith varying boundary
condition): /? = 1 ; k = 1 .....................................................................................
61
Steady-state curves (constant complex perm ittivity cases): (3 = 1; k = 1;
Tbw = - . 5 ................................................................................................................
61
4.6
M elting front position S vs. derivative S': /3 = 1; k = 1.................................
66
4.7
Steady-state curves (constant and variable complex perm ittivity cases):
(3 = 1; k = 1; Tbw = - . 5 ......................................................................................
67
M elting front position S vs. tim e t (constant complex dielectric perm it­
tivity): {3 = 1; k = 1.............................................................................................
74
Melting front position S vs. tim e t (tem perature dependent complex
dielectric perm ittivity): P=400; f3 = 1; k = 1 ...............................................
74
4.4
4.5
4.8
4.9
4.10 M elting front position S vs. tim e r : P = 600; (3 = 1; k =
1;
= —.5. .
76
4.11 M elting front position S vs. tim e t: P = 600; /? = 1; k = 1 ; Tbw = —.5. .
76
x
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 1
IN T R O D U C T IO N
1.1
H ydrodynam ic S tab ility
W hen most fluids are heated, their densities decrease. If th e tem perature profile is
such th a t a more dense fluid lies over a less dense fluid, then we have a poten­
tially unstable situation in the presence of a gravitational field.
In the present
study we wish to determine the effects of microwaves on the tem perature profile of a
given dielectric fluid and, consequently, their effect on the velocity field of the fluid.
This problem arises in separating emulsions Fang and Lai (1995), and in coating
m aterials. In the latter case, powders are distributed on a surface and are heated
with microwaves until they m elt and smoothly coat the m aterial. These applications
m otivate the use of both rigid-rigid and rigid-free boundary conditions. (Note th a t
we do not consider interfacial deformations in the present study.) One of the main
param eters of concern is th a t connected to incident power levels necessary to induce
convective m otion within the fluid layer.
Much work has been done in th e area of convection driven by internal heat
sources.
The problem arises in geophysics and may be of im portance in meteo­
rology. Pellow and Southwell (1940) have studied the hydrodynam ic stability of a
fluid subject to a linear tem perature profile. Sparrow, Goldstein and Johnson (1964)
studied the effects of imposing a nonlinear tem perature profile, due to a uniform
volum etric heat source, on th e stability of a fluid layer. T he theory treats convective,
insulated and fixed boundary conditions for th e tem perature. Roberts ([10],[11]),
studied th e convection p atterns in a horizontal layer of fluid driven by a constant
internal heat source. Using energy m ethods he constructs an iterative scheme to
study the nonlinear stability characteristics of th e fluid and resolves degeneracies left
by linear stability theory. W atson (1968) studied th e effects of introducing a uniform
volum etric heat sink into an infinite horizontal fluid layer on the stability of th e
1
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2
fluid w ith fixed boundary conditions for the tem perature. She further investigates
the possibility of formation of double layers of convective cells for large values of
the rate of heat loss by the stable sublayer of fluid at the top and concludes that
such a formation does not occur if the tem perature a t the upper surface is less than
th a t a t the lower surface. Studies were performed for the rigid-rigid and free-free
boundary cases. Yucel and Bayazitoglu (1979) studied the effects of a nonuniform
volumetric heat source on the stability of an infinite horizontal fluid layer. The source
is induced by absorption of incident radiative energy penetrating into the fluid. They
impose convective boundary conditions on the tem perature and rigid-rigid boundary
conditions on the fluid.
In studies focusing on hydrodynamic instabilities within the bulk fluid, the
source in the heat equation was chosen either for m athem atical convenience or as
a simple model of radiative heat transfer. There are many applications in which
a microwave source becomes the heat source for fluids such as w ater and ethanol
and solids such as ceramics.
This is owed to the fact th a t for these and other
m aterials, microwaves can be an effective means for depositing energy rapidly and
selectively within the m aterial. The mechanisms for this deposition are ohmic and
dipolar heating (see Von Hippel (1954), Ramo (1984), Feynman (1964)). The use
of microwaves in the processing of materials, ie. heating, sintering, melting, is an
active area of research. A study of microwave heating a cylindrical vessel of fluid
w ith a free surface has been carried out by Stuerga and Lallemant (1993a,b), (1994),
and Stuerga, Lallemant and Steichen-Sanfeld (1994). Experiments are presented
to investigate microwave heating of a confined fluid under reduced pressure (free
surface) and make conclusions about boiling, superheating and explosions occurring
in a fluid heated with microwaves. Their work further suggests th e ability to control
the shape of the spatial therm al profile within the subject fluid through microwave
irradiation. Stuerga, Lallemant and Zahreddine (1993) studied th e interaction of
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3
microwaves with water. This study gives evidence th a t therm al runaway can occur
in w ater at appropriate microwave power levels. This is owed to the tem perature
dependence of the complex perm ittivity of water which can give rise to resonance
of the electric field within this media and, hence, can give rise to bistable solutions
for steady-state tem peratures.
M. Nachman and G. Turgeon (1984) studied the
interaction of microwaves in multilayered material. Specifically, they studied the
heating patterns th at evolve within a three layered m aterial and investigated the
ability to preferentially heat a specific layer by introducing a flat reflector to create
standing waves within th e media.
In chapter 2 we consider in detail th e microwave interaction with a fluid and
the induced source term.
Specifically, we see th a t changing the impinging wave
frequency or the slab thickness alters th a t position where the basic tem perature
profile attain s its maximum value within the fluid medium. We conduct a study of
the hydrodynamic stability of a fluid layer heated with a microwave source. We make
the assumption th a t the dielectric attrib u tes of the subject m aterial are independent
of tem perature. The unperturbed solutions are found and the linear stability of these
basic solutions is investigated. The analysis suggests minimum power requirements
necessary to generate convection w ithin th e fluid. We also find th a t this minimum
power undulates as the thickness of the fluid slab or frequency of th e microwave is
altered.
In chapter 3 we develop a model which incorporates the tem perature dependence
of the complex perm ittivity of a m aterial. T he unperturbed solutions are found.
We specifically focus on the steady-state tem perature distributions. As will be seen,
a therm al runaway phenomena can occur in certain fluids. We study th e hydrodynam ic stability of a fluid layer heated w ith a microwave source. It will be shown
th a t the power levels needed to achieve an onset of convection within th e fluid layer
fall far below those necessary to cause therm al runaway within the fluid.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
4
1.2
A pplications o f M icrowaves in th e M eltin g o f M aterials
As m entioned in the first p art of our introduction, microwaves have potential
applications involving the melting of solid m aterials such as in welding processes,
rem oval/application of materials to /fro m surfaces and thawing frozen foods.
In studying such processes it is im portant to develop an understanding of the
mechanisms describing the interaction of the microwave source w ith m aterials, as well
as an understanding o f the mechanisms involved in phase-change front form ation and
propagation. In applications such as the coating application mentioned previously
a change of phase is involved in melting the solid to coat a surface. It is desirable
to understand this m elting process in order th a t it can be carried out in an efficient
m anner, (i.e. minim al expenditure of power is involved or a faster turnaround time
is achieved). The use of microwaves in melting, may be useful in such applications
as coating application/rem oval in th a t they provide a more uniform heating of many
m aterials th an do conventional heating sources. For m any m aterials, microwaves may
significantly reduce th e likelihood of surface dam age in the coating removal process.
For example, in the removal of paint from wooden surfaces, the surface is less likely
to be scorched if microwaves are used than if a conventional heating gun is used.
Many studies b o th theoretical and experim ental have been made in the area
of phase front propagation (see B. A. Boley (1961), I. C rank (1957), S. J Citron
(1960) and others). O ne underlying assum ption in m any of these works is the latent
heat balance (Stefan) condition along with a Dirichlet condition on the tem perature
im posed a t the m elting front. For a more comprehensive investigation of Stefan
problem s see Hill and Dewynne (1987), C rank (1984) and Rubinstein (1971). In
short, the Stefan condition states th a t the m elting front velocity is directly propor­
tional to th e difference in the energy fluxes on b oth sides of th e m elting front. In many
of these studies the source/sink is imposed through th e boundary. O ther works by
C rank (1984) deal w ith moving fronts induced by volum etric sources using enthalpy,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
5
front fixing and other methods. Goodman (1958) approxim ated the propagation of
a melting front by assuming a polynomial representation of the tem perature in the
liquid region and using a heat balance integral to construct a governing ordinary
differential equation for th e melting front. Huang and W ang (1994) investigated the
propagation of a freezing front during m etal casting w ith mold. They employ inverse
m ethods to deduce a suitable boundary condition a t the m old/casting interface and,
with this estim ated boundary condition, im plem ent the enthalpy m ethod to deduce
the trajectory of th e melting front. Sethian and S train (1991) present an efficient
algorithm for predicting crystal growth and dendritic solidification in supercooled
media.
In chapter 4 we investigate the melting of dielectric solids using microwaves. We
will trea t the one-dimensional case. Specifically, we investigate the trajectory of the
melting front in tim e and raise issues regarding feedback control on the power source.
The model will incorporate Maxwell’s equations along w ith an energy equation for
both the liquid and solid regions. As we will see, there is a strong coupling between
the melting front position, the electric field strengths of the liquid and solid regions
and the tem perature field of the liquid region. The Stefan condition is employed
to model the movement of the melting front. Having posed the full m athem atical
model, th e dynamics and steady states of the system are investigated.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CH APTER 2
H Y D R O D Y N A M IC STABILITY OF A FLU ID LAYER: C O N ST A N T
DIELECTR IC P E R M IT T IV IT Y CASE
2.1
T he M odel
We consider a viscous incompressible fluid of constant therm al conductivity k t
and viscosity p., bounded between two infinite parallel plates a distance d apart.
The initial tem perature of the fluid, 0O, is assumed to be the same as th a t of the
surrounding space, and the corresponding initial fluid density is denoted by p0.
For simplicity, the bounding plates are taken to be transparent to th e im pinging
microwave field. The geometry is shown schematically in Figure 1.
The equations governing this microwave heating problem are th e tim e dependent
Maxwell’s equations for a moving medium, the energy equation, and the NavierStokes equations. In our model we assume th a t the fluid velocities generated are
much smaller than the speed of light. Then the relativistic correction term s to the
electromagnetic equations are negligible and Maxwell’s equations reduce to their
standard form in a stationary medium.
The second assum ption we make is the
charge neutrality of the fluid. Thus, the Lorentz force in the m om entum equation
for the fluid is absent. Finally, we assume th a t the time required for the fluid or
heat to diffuse a wavelength is much larger than a microwave period. T he later
is O (10-10) seconds for commercial generators. This last assum ption allows us to
average all the governing equations over a microwave period (see appendix A). The
resulting equations are a time harmonic version of Maxwell’s equations and the
tim e dependent Navier-Stokes and heat equations. T he later contains the averaged
microwave source term .
W ithin this framework, we assume th a t a plane time-harm onic electrom agnetic
wave of frequency ui impinges normally upon the fluid layer which fills th e region
0 < z' < d. A portion of this wave scatters from the interface z ' = d, a portion
6
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
L = electromagnetic wave length
F ig u r e 2.1 The geometry of the hydrodynamic problem.
penetrates the layer and heats the fluid, and the remaining portion is transm itted
through th e interface z' = 0. In the free space regions z' > d and z' < 0, th e electric
field is given by the real p arts of
( 2 . 1)
E = E o T e -i(-k'z'+ut,)XL,
z' < 0,
(2 .2 )
respectively, where Eo is th e strength of the incident field, k' = uj/ c, c is th e speed of
light in free space, T is the transmission coefficient, and
Both T and
7
7
is th e reflection coefficient.
are to be determined.
The electric field which penetrates th e fluid and interacts with it is given by
the real p art of E = [E'2{z')exp(—zwf')]]x, where E 2 satisfies
(2.3)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
In this equation k[ = (w/c) y/e^, e0 is the perm ittivity of free space, ex is the perm it­
tivity of th e fluid and er = ei/e 0 is the relative perm ittivity. If the fluid is purely
dipolar, then a0 is the im aginary p art of the dielectric constant. On the other hand,
if th e fluid is dipolar with conductive losses, then <7o is the effective conductivity,
Kriegsmann (1993). In either case £ = ^
is a measure of how much power from
the incoming wave is absorbed by the fluid; see Hippel (1954), Ramo, W hinnery and
Van Duzer (1984). Implicit in the definition of k[ is our assum ption th a t the fluid is
non-magnetic. We also assume here th a t all the electrical param eters, ju st defined,
are independent of tem perature.
From the continuity of the tangential electric and m agnetic fields a t z' = cl
and z' = 0, we deduce th at E 2 and its derivative are continuous there Ram o et al.
(1984) (see section 3.14), Kriegsmann (1992). Combining this fact with (2.1-2.2) and
elim inating T and
7
, we find th a t E 2 satisfies the boundary conditions
+ ik'E'2 = 0,
z' = 0,
z' = d,
(2.4)
(2.5)
Equations (2.3-2.5) fully determ ine the electric field w ithin the fluid and outside the
fluid.
The electric field described above acts as a source of energy th a t heats the
fluid. The fluid equations to be solved consist of th e Navier-Stokes equations for
heat conducting fluids coupled to a tem perature equation which includes an energy
supply due to th e electric field. Denote dimensional fluid velocities, tem perature and
pressure by u ' = (u\,u'2,u'3), 9' and P ' respectively. M aking the usual Boussinesq
approxim ation, p = p o ( l—
~ $o)) where a is the coefficient of therm al expansion
of th e fluid, (see Drazin & Reid (1981)) leads to the following system:
_
Po'~Dt = ~ 9Po{1 ~ a{6' ~ 9o))6i3 ~ V P ' +
V • u ' = 0,
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
( 2 .6 )
(2.7)
where x ' = (x ', y', z') and A = ^
+ -§ji- The boundary conditions are those
of no slip for fluid velocities a t the solid walls
u[ = u'2 = U3 = 0,
z = 0,d,
(2.9)
In the case of a free upper surface th e boundary condition a t z = d becomes
du\
du'o
du'o
du'o
a7+ a?=a? + #
.
,
= 0' “’ = 0'
^
fn
(2'10)
and convective conditions for th e tem perature a t the walls
BB*
k l a i = h{9' ~ °o)’
z = 0>
(2’n )
BB*
k t =
~ h{-e' ~ 0o)’
z=
(2’12)
In equations (2.6-2.12) above g is th e acceleration due to gravity, Cp is the specific
heat of th e fluid and h is the h eat transfer coefficient m easuring loss of heat due to
convective cooling. The source term in (2.8) is due to the microwave interaction with
the fluid.
The governing equations are made dimensionless by scaling distances with the
slab thickness d, time by the diffusive tim e scale d2/ k where k = k t/{poCp) and the
electric field by the incident am plitude E q. New variables are introduced as
x = x '/d ,
t = Kt'/(d2),
P = d2/(K 2p0)P',
u = du'/«
9 = (-1 + £ ) / *
(2.13)
New electric field variables are introduced as
E = E ' / ( e - ik'E 0),
k = k'd,
k x = k[d,
(2.14)
Substitution of (2.13),(2.14) into (2.3-2.5) and (2.6-2.8 ) gives th e dimensionless
system given below. T he nondimensional groups which enter into th e dynamics
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
10
are
the Rayleigh number, R, the P ran d tl number, P r,
the Biot number,
a nondimensional measure of the incident microwavepower
0, and
denoted by x- These
groups are given by:
R = {o c g d 390)/{Ku),
P t = u/ k ,
0 = h d / k t,
X = a 0\E0\2d2/{ 2 k t90),
(2.15)
These scalings produce the nondimensional system:
^
j2 p
— + k 12(l + i O E = 0,
(2.16)
V • u = 0,
(2.17)
= - V ( P + gzd3/ k 2) - R x P r f fk + P r A u ,
(2.18)
^ - = A 9 + \E\2,
(2.19)
dE
—— I- i k E — 0,
2 = 0,
dz
dE
—
ik E = —2ik,
2 = 1,
dz
B9
— = 09,
2 = 0,
= -0 9 ,
^
Free upper surface
ay
ox
where k is the unit vector in the
2
ox
(2.21)
(2.22)
z = 1,
Ui = u2 = u3 = 0,
Rigid —rigid
(2.20)
oy
(2.23)
2
= 0, u3 = 0,
= 0,1,
2
= 1,
(2.24)
(2.25)
direction. It should be noted th a t in th e basic
state (ie. u = 0), the governing equations for th e tem perature 9 and the electric field
E are independent of the flow and are not affected by the Boussinesq approxim ation.
2.2 B asic States
T he base flow velocities are u = 0 and equation (2.18) gives th e corresponding
pressure distribution. Since equations (2.16) and (2.19) are decoupled, as long as
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
11
f is independent of tem perature, we can explicitly determ ine the electric field and
hence the tem perature. The solution for the electric field is given by:
2ik[T c o s (r(l - z)) - ik s in ( r ( l - z))\
sin(r)(r2 + k2) + 2 ik r cos(r)
’
-
[
’
where T = k\-^(l + i£). The corresponding tem perature field is obtained by substi­
tuting (2.26) into (2.19), integrating with respect to z, and applying 2.22 and 2.23;
the result is
0(z) = - J qZ\E\2(z - s)ds + {fiz + 1) j ^ r p ) f Q \ E \ \ l + P(1 ~ s))ds,
(2.27)
Note th a t the electric field (2.26) depends on k (this param eter is the ratio of the slab
thickness to the microwave length) and the param eter £ through the param eter T.
In Figures 2.2a-d we give typical undisturbed tem perature distributions within the
slab for different microwave frequencies k for fixed T and unit Biot number, ft — 1.
(We mention th a t these tem perature profiles are those needed to achieve an onset of
instability within fluid, as is explained in Section 3.) As can be seen from th e profiles
in figure 2 .2 , the position of maximum tem perature within the slab is a complicated
function of k - in fact it is found to oscillate with k as we will see later.
2.2.1
A Lim iting Case k -> 0 (Low Frequency)
The case of long microwaves, k «
1, can be analyzed asymptotically. Assuming
th a t er is of order one we can expand in powers of k where ki = erk to approxim ate
the electric field as E = E„ + E \k + E2k2
S ubstituting this into (2.16) and
utilizing th e boundary conditions ( 2 .2 0 ),( 2 .2 1 ) yields the following system:
C?2 (jgo+^ f
dz
1
+ " ' ) + erk 2{1 + iO {E 0 + k E x + •••) = 0,
+ "
- ik (E 0 + k E l + ---) = —2ik,
d(Eo+A E ' + • • •) + i k {Eo + k E l + . . . ) = 0,
2
(2.28)
z = 1,
(2.29)
= 0,
(2.30)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
12
Fig. b (k= 0.705)
Fig. a (k= 0.205)
0.8
0.8
0.2
0.2
1.7
1.007
1.8
1.9
dimensionless temperature
Fig. c(k= 1.205)
1.008
1.009
dimensionless temperature
Fig.d(k= 1.705)
1.01
0.8
0.8
0.6
0.4
0.2
0.2
1.014
1.018
1.016
dimensionless temperature
1.02
1.0055
1.0045
1.005
dimensionless temperature
1.004
F ig u r e 2 .2 Typical undisturbed tem perature profiles; the physical param eters are
0 = 1 ; er = 2.54 ; £ = .25.
The leading order system is:
d2E 0
= 0,
dz2
dE0
= 0, z =
dz
dE0
■= o, z =
dz
(2.31)
1,
(2.32)
0,
(2.33)
Solving th is system we find E 0 = B where B is any constant. To solve for B we
m ust look to the next order system:
<f E l
=
dz2
0
(2.34)
,
dE,
- —ikEp = —2 ik, z = 1,
dz
dEl
+ ik E 0 = 0, z = 0,
dz
Solving for E x, we find E i = C z + D so ^
(2.35)
(2.36)
= C. We now utilize th e boundary
conditions to construct an algebraic system th a t lets us solve for B :
C —ik B = —2ik,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.37)
13
C + ik B = 0,
(2.38)
Hence B = E a = 1 so E = 1 to leading order. Substituting this expression for E
into (2.27), we obtain the following:
0(z) = - z 2/ 2 + z / 2 + 1/(2/?),
This is the familiar parabolic tem perature distribution th a t would
(2.39)
be obtained from
constant volum etric source heating (see Introduction).
2.3
Linear Stability
The linear stability of the undisturbed states described above, can be analyzed by
adding small perturbations u ', 9', P ' and E' and linearizing (2.16-2.25) w ith respect
to prim ed quantities to obtain:
V • u ' = 0,
(2.40)
dn*
^ = - V P ' - R x P r d 'k + P r A u ',
(2.41)
% +
at
dz
= AS".
(2.42)
Taking the curl of (2.41) twice, we obtain the following equation involving th e vertical
velocity com ponent of the fluid and the tem perature:
(/I
where A !
=
= R x P r A iO ' + P r A V ,
+ Jp-. This equation along with (2.42) willbe used
(2.43)
to determ ine
the stability of the system. Note th a t w ith th e eigenfunctions w' and O' known the
continuity equation (2.40) along with (2.41) yields the velocity com ponents v! and
v' as well as th e perturbation pressure P ' .
E quations (2.42) and (2.43) comprise a coupled system of PD Es which can be
analyzed by introducing the norm al modes
O' = Q ( z ) f( x , y)eat,
w = W ( z ) f ( x , y)ea\
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.44)
14
into (2.42) and (2.43). For example, the energy equation (2.42) becomes
r
, ezW
0
D 2Q 1
A ,/
0 J
/ ’
which on introduction of a separation constant a2, yields the Helmholtz equation
A i/ + a2f =
0
for the horizontal dependence of th e tem perature eigenfunction.
Proceeding as above and using the momentum equation (2.43) also, yields the
following eigenvalue problem depending on the wavenumber a which arises from the
Helmholtz equation and which measures the wavelength of linear perturbations:
{D2 - a2)(P r (D 2 - a2) - a ) W - a2R x P r Q = 0,
(2.45)
(D 2 - a 2 - a ) e - W D 0 = 0,
(2.46)
W = D W = 0,
Rigid - rigid
2
W = D 2W = 0,
Rigid - free
= j3Q,
^
~ = -/?© ,
dz
= 0,1,
(2.47)
z = 1,
(2.48)
z = 0,
^ = 1,
(2.49)
(2.50)
M arginal stability of the above system is given by a = 0. Although exchange
of stabilities has not been proven in th e present problem (the difficulty arises due to
the fact th a t
is non-constant in th e slab), it is strongly suggested numerically.
2.4
N um erical Solutions
Due to th e complicated form of the background tem perature variations (see (2.26)
and (2.27)) a closed form solution of th e eigenvalue problem has n ot proven possible,
and so we proceed numerically.
later.
However, some lim iting cases are considered
The problem was solved numerically using a M atlab O D E solver (ode45)
in conjunction with a shooting m ethod. Before proceeding with a description of our
numerical results, we make some comparisons w ith o ther related investigations.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
-2
-4
-6
R'X = 104430
— R'X = 125000
R*X = 150000
-8.
wave # (a)
F ig u re 2.3 Growth rate curves, rigid-rigid case; (5 = 1; k = 1; er = 2.54; P r = 17;
£ = .25.
Roberts (1967), has considered convection problem s w ith uniform internal heat
sources (i.e. producing parabolic undisturbed tem perature profiles), with isothermal
boundary conditions on the upper face and an insulating condition on the lower face
- these conditions translate, in our notation, to 6 =
0
at z =
1,
and jj| =
0
at
2
=
0
and similar expressions for the perturbation. This scenario can be recovered from
our system by (i) letting the Biot number /? —* 0 0 on th e upper face (this gives 0 = 0
at
2
= 1 to leading order),(ii) letting (5 —> 0 a t
face), and, (iii) letting k -¥
0
2
= 0 (giving ^ = 0 a t the lower
to obtain th e background tem perature distribution
given by (2.39). O ur code reproduces the neutral stability results of R oberts noting
th a t his Rayleigh number, in our notation, is given by R xAn additional test of our numerical procedure is possible in the lim iting case
of /3 —>• 0 on b oth the upper and lower faces, along w ith k —» 0 to produce uniform
volumetric heating. This is a special case addressed in th e stability study of Yucel
and Bayazitoglu where th e attenuation coefficient is zero. O ur results are in complete
agreement with theirs. In particular we obtain the same critical value of th e Rayleigh
number R x = 37328.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
16
2.4 .1
A L im itin g C ase (Low W av e N u m b e r P e r tu r b a t iv e M o d e s )
It is seen from Figure 2.3 th a t in the long wave lim it a —» 0, the flow is stable with
a decay rate independent of R x ■ The numerical value found for the set of physical
param eters of Figure 2.3, is —1.707. This result can be confirmed analytically by
an asym ptotic expansion for small a. In the lim it a —> 0, th e system (2.45-2.50)
becomes, to leading order,
{D4 - - ^ D 2) W = 0 ,
(2.51)
(D2 - a)Q - DOW = 0,
(2.52)
W = D W = 0,
z = 0,1
(2.53)
DQ = 0 0 ,
z - 0,
(2.54)
DQ = - 0 0 ,
z =
(2.55)
1
M ultiplying (2.51) by W* and (2.52) by ©* where stars denote the complex
conjugates, and integrating by parts we arrive at the following equations:
[ l \DQ\2d z - a [ l \Q\2d z + [ ' D(9)Q*W dz = 0,
Jo
Jo
Jo
(2.57)
First we use (2.56) and (2.57) to show th a t a must be real. If W is non-zero,
then (2.56) is sufficient to determ ine th a t a is real. If W (z ) = 0, a is seen to be real
from (2.57). In addition th e expressions (2.56) and (2.57) show th a t if a = 0 then
b oth W ( z ) and 0 (z ) are zero resulting in a trivial solution.
It has been found numerically th a t the W eigenfunction tends to zero as a
decreases, while the corresponding 0 eigenfunction is non-trivial. This behavior is
illustrated in Figures 2.4a,b where W {z) and 0 (z ) are plotted for different decreasing
wavenumbers and R x = 125000, the other param eters being those of Figure 2.3.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
17
Fig. 4a
600
400
sr 200
-200
0.1
2.5
02
0.3
0.4
0.5
0.6
height (2 )
Fig. 4b
0.7
0.8
wave # (a) -1
- • wave # (a) ».5
-
wave 4 (a) «.1
1.5
0.1
0.2
0.4
0.7
0.9
Figure 2.4 Eigenfunctions W (z ) (top) and 0 (z ) (bottom ); /? =
P r = 17; £ = .25.
1;
k = 1; eT = 2.54;
W riting W ( z ) = 0 + . . . , then, the system (2.51-2.55) is decoupled and can be solved
exactly for
0
to obtain the following eigen relation for a:
t a n ( \ /—cx)(/32/ \ / — — \ [ —&) + 2/3 = 0,
(2.58)
The eigenrelation (2.58) was solved numerically to determ ine the largest eigenvalue
for a to be -1.70, in complete agreement w ith the numerical solutions presented in
Figure 2.3. Note th a t this value is independent of R x as is dem onstrated in Figure
2.3.
2.4.2
A L im iting Case (S tability o f M odes for Large C ritical Power
Levels)
As can be seen from Figure 2.6, an unstable mode exists if the (modified) Rayleigh
num ber is large enough. In fact, asym ptotic m ethods can be used to analyze the most
unstable mode as R x —>oo as presented next. Consider equations (2.45)-(2.46) with
R x asym ptotically large and the wavenumber a of order one. The leading order
balance leads to an inviscid system away from the walls. More formally, a balance
of term s in (2.45)-(2.46) leads to the expansions
W = (R x )1/2W 0 + . . . ,
0 = 0o + . . . ,
a = (R X) 1/2ct0,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.59)
18
where cro(a) is the scaled eigenvalue to be found.
These expansions lead to the
inviscid leading order eigenvalue problem (after elim ination of 0o between (2.45)
and (2.46)):
a 20 {D2 - a2)W 0 - a2f - W 0 = 0,
az
Wq=
0
,
on
z=
0, 1,
(2.60)
We note th a t in the case of a parabolic basic tem perature profile 0, the system (2.60)
is analogous to the stability of Couette flow in the gap between concentric cylinders
at large Taylor numbers (see Drazin & Reid
(1981)), which can be solved in term s
of Airy functions; in th a t special case it is
found th a t oois bounded above by a
constant, Co say, and cr0 —>■c0 as a —> oo. The m ost unstable mode is for short waves
and viscous boundary layers are required to reduce this to zero as th e right hand
neutral branch of the stability curves is approached (see Figure 2.6, for instance).
In the case of general wavelength microwaves the eigenvalues m ust be com puted
numerically for different values of a. Similar trends as in the parabolic profile case
are found. Typical results of solutions of (2.60) are shown in Figure 2.5; th e eigen­
function W q( z ) is shown for wavenumbers a = 1,100,1000 for the physical param eters
indicated. It is seen th a t as a increases the mode is concentrated near z = 1 and is
zero away from it. This can be quantified by constructing W KB solutions of (2.60) as
a —»• oo; a turning point exists in these solutions near z = 1 and solutions valid below
and above the tu rn in g point can be constructed by standard m atching m ethods. The
m atching provides an upper bound on th e leading order eigenvalue cr0,found to be
<7 o —>• - m in |
,
0
For th e param eters given in Figure 2.5, the
< z <
1.
result is a \ —>
.1277 as a —> oo. See
A ppendix C.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
19
z
400
•*■100; slgm«0B.O42
4
xto”
z
|2
z
F ig u r e 2.5 Asym ptotic eigenfunctions a t large Rx', P = 1; k =
2.5
1;
er = 2.54; £ = .25.
R e s u lts
Stability characteristics were obtained by fixing th e physical param eters including
the microwave power (in non-dimensional term s this is x) which in tu rn fixes the
power Rayleigh number R x ■ For different wavenumbers a, th e eigenvalue a, which
is complex in general, was computed by iteration. O ur numerical results produce
real values of a alone, indicating th a t there is no propagation of disturbances down
the slab. In addition these results suggest th a t neutral curves can be traced by
specifying a = 0. A typical set of results is given in Figure 2.3 for different power
Rayleigh numbers and physical param eters (shown on the Figure) corresponding
to ethyl alcohol. As can be seen, larger values of R x support a band of unstable
modes which narrows as the modified Rayleigh num ber decreases - see the curves
corresponding to R x = 150000 and 125000 respectively. As th e modified Rayleigh
num ber is decreased further, it reaches a critical value below which th e flow is linearly
stable to all wavelengths. The critical value for th e set of param eters of Figure 2.3
is found to be 104430 and the corresponding variation of a w ith wavenumber is also
given. It can be concluded, therefore, th a t the qualitative effect o f increasing the
microwave power is an enhancement of the instability.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
20
Next we use typical stability plots, as in Figure 2.3, to present neutral curves
which have a = 0. For a fixed R x it is seen th a t there are a t m ost two values
of the wavenumber where the flow is neutrally stable, and below a critical value of
R x , denoted by R x c, the flow becomes linearly stable to all wave numbers. N eutral
curves of R x against a for different Biot numbers (other param eters being fixed) are
considered first. In Figure 2.6 we show neutral curves for /3 = 0.01 and 1.0 when
the non-dimensional param eter k = 1.0. The unstable and stable regions are located
above and below the given neutral curve respectively. In fact each neutral curve
provides a critical minimum point (ac, Rxc)\ flows corresponding to Rayleigh numbers
below R xc are stable. The variation of (ac, Rxc) in param eter studies provides a lot
of inform ation about the stability characteristics and is considered later. Physically,
(3 may be associated with the rate at which heat is convected away from the slab and
hence an increase in the Biot number has a stabilizing effect on th e system. T h a t is,
for a given value of R x the window of unstable modes contributing to the disturbance
in the fluid decreases. In the case of a free upper surface the m arginal curve drops
indicating th a t the fluid is more prone to become stable. These observations indicate
th a t the more constrained a system, the greater its stability.
One m ain difference between background heating profiles produced by microwave
radiation w ith those produced by constant volumetric heat sources, is th a t in the
former case the height (z) where the background tem perature is a maximum varies
w ith k, whereas it is fixed a t z = 1/2 when k tends to zero (see (2.39)). This has
an im portant influence on the convection p atterns th a t emerge since the point of
m aximum tem perature separates regions of negative therm al gradients above and
positive ones below. It can be expected, then, th a t in regions where th e therm al
gradient is positive the flow is stable and the resulting perturbation velocities are
small.
Convection patterns in the form of rolls should have centers which are
displaced vertically upwards.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
21
2.S
x 10
— rigid-rigid'
- • rigid-rigid'
Biot # « .01
Biot f * 1
rigid-free ca se Biol • ■ 1
1
W ave * (a)
F ig u re 2.6 N eutral stability curves in the
X
10
plane; k = 1; er = 2.54; ; f = .25.
X 10
05
1
1.5
2
2.5
ratio of slab thickness to microwave length (k)
3.5
0.5
1
15
2
25
ratio of slab thickness to microwave length (k)
3.5
*33
.£
&
£
‘8 1
0
1■3
19
1.5
2
2.5
1
ratio of stab thickness to microwave length(k)
3
3.5
F ig u r e 2 .7 N eutral curves rigid-rigid case: (3 = 1; er = 2.54; £ = .25.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
22
11-------------------i------------------ 1-------------------- 1-------------------1------------------- 1---------------------1------------------1
0
•5
X
0.5
1
1.5
2
2.5
ratio of slab thickness to microwave length (k)
3
3.5
I_ _ _ _ _ _ _ _ i _ _ _ _ _ _ _ _ _ _ _ _ i _ _ _ _ _ _ _ _ _ i _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ i_ _ _ _ _ _ _ _ _ _ _ _ i_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ i _ _ _ _ _ _ _ _ _ _ p
0
0.5
1
13
2
23
ratio of stab thickness to microwave length(k)
3
3.5
F ig u re 2.8 As in Figure 2.7 but rigid-free boundary conditions.
We now quantify these observations by presenting numerical solutions of the
eigenvalue problem.
Figure 2.7a gives the variation of the critical values of R x c
w ith k, the ratio of slab thickness to microwave length. Figure 2.7b is a plot of the
critical wavenumber ac versus k, and Figure 2.7c presents the variation of th e height
where th e background tem perature achieves its maximum, zmax say, with k. The
critical values for R \ are those values a t which the onset of instability is achieved
and exactly one perturbation mode corresponding to the critical wave num ber ac
becomes excited. It can be seen from these results th a t th e variations with k of the
three quantities plotted in Figure 2.7a-c are roughly in phase. This can be explained
physically as follows: Consider Figure 2.7c first. A decrease in the height zmax implies
th a t there is a larger unstable layer of fluid residing near th e upper boundary. The
lower stable region is almost stagnant (see later also). The size of convection cells is
expected to scale with the size of the unstable layer, and so decreasing zmax increases
the size of the cells which are proportional to 1 /ac. This argum ent indicates th a t the
variations of ac and zmax w ith k should be roughly in phase, and this is supported by
the numerical results. Figure 2.7a also indicates an in-phase behavior between R x c
and zmax■ This can be attrib u ted to the fact th a t a decrease in zmax provides a larger
unstable body of fluid to be convected, which could be achieved with a lower power
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
23
input, i.e. a lower Rxc- Again, the num erical results produce in-phase variations
w ith th e exception of small values of k. We note th a t th e non-monotonic behavior in
Figure 2.7a, for example, is capable of producing th e same critical Rayleigh num ber
for three different microwave lengths (see Figure). This mechanism can be used to
obtain sim ilar convection cell sizes for different values of k.
Figures 2.8a-c show
sim ilar plots for a lower rigid and free upper boundary. T he critical values of R x
are lower as might be expected since tangential slip is now perm itted a t the upper
slab face and the system is less constrained. Finally, we observe th a t in both cases
the overall trend of R x c with k is increasing. Since R x scales w ith E %, the power of
the incident field, and k scales w ith th e slab thickness d, this trend indicates th a t
thicker fluid slabs require more power to become unstable. The local minimum at
k ~ 2 occurs because the m agnitude of th e electric field deduced from Figure 2.7
has a m axim um there.
This is a resonance phenomenon. Thus, th e strength of
the incident field and the value of R x c required for the onset of an instability are
reduced. Similarly, the local maximum a t A: ~ 3 occurs because the m agnitude of
the electric field is minimum. This is an anti-resonant behavior. Thus, the strength
of the incident field and the value of R x c required for the onset of an instability are
increased.
In closing this section we present a typical flow which shows convection cells at
the onset of instability. Two values of k (.705, 1.355) are chosen, yielding a zmax less
th a n 1 /2 and greater than 1/2 respectively. T he Biot num ber in both cases is one,
and the m aterial param eters are those for ethyl alcohol. C ritical wavenumbers and
Rayleigh numbers are obtained numerically as described earlier, and th e convective
cells shown in Figures 2.9-2.10 and were generated by choosing f ( x , y ) = cos (ax)
and utilizing the com puted eigenfunctions. Also shown on th e figures are the corre­
sponding background tem perature profiles (note th a t these have been scaled up for
illustrative reasons). The figures support our earlier observations of Figures 2.7a-c,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
24
1 .2
Pertu rb atio n Flow R eid
r
-0-0.2.5^--------01------0.51------- 11--------1.51--------21------2.51---------31------'
3.5
Horizontal position
F ig u r e 2.9 P erturbation flow field a t onset; base tem perature profile x 10 super­
imposed on flow field. (3 = 1; k — .705; er = 2.54; f = .25; critical wave num ber ac
=2.45.
P erturbation Flow R e id
1.2
0.8
0.6
S
0.4
0.2
0.5
1.5
2.5
3.5
Horizontal position
F ig u r e 2.10 Perturbation flow field at onset; base tem p eratu re profile x 10 super­
im posed on flow field. 0 = 1; k = 1.355; er = 2.54; f = .25; critical wave num ber ac
=4.401.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
25
Perturbation Flow Field
1.2r
-Qo '
-0.5
1------ 1---------1------- 1
0
0.5
1
----------------------
1.5
1------
2
*----------------- 1----------------------- '
2.5
3
3.5
H orizontal position
F ig u r e 2.11 Perturbation flow field at onset; base tem perature profile x 10 super­
imposed on flow field. /? = 1; k = .71; er = 2.54; £ = .25; critical wave number ac
=1.99.
Pertu rb atio n Flow R eid
0.8
0.6
5 0.4
0.2
0.5
1.5
2.5
3.5
H orizontal position
F ig u r e 2.12 P erturbation flow field at onset; base tem perature profile x 10 super­
imposed on flow field. (3 = 1; k = 1.36; er = 2.54; f = .25; critical wave number ac
=3.22.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
26
and it is seen th a t the cell-size in Figure 2.9 is larger th an th a t of Figure 2.10; in
addition the cell centers lie above zmax and flow velocities below zmax (i.e. in the
stable region) are smaller than those in the unstable region. Figures 2.11 and 2.12
show representative flows for the rigid-free boundary case. These graphs support the
results of figures 2 .8 a-c.
2 .6
D iscu ssio n
We have considered convective instabilities in fluids heated by microwaves.
For
constant electrical conductivities it has been found th a t the induced background
tem perature gradients produce convective cells when a critical Rayleigh num ber is
exceeded. Microwave heating is equivalent to volumetric heating but th e profiles
produced are more complicated since they result from th e power supplied by the
non-uniform electric field within the slab.
One common feature, however, is the
presence of a tem perature maximum a t some height in the slab and consequently a
region of stable fluid (cooled from below) and a region of unstable fluid (heated from
below). An interesting feature of the microwave produced background tem peratures,
is th a t th e position of maximum tem perature can be varied by varying the microwave
frequency; the variations are oscillatory with the non-dimensional param eter k (see
earlier definition).
It has been found numerically th a t critical Rayleigh numbers
as well as critical wavenumbers are qualitatively in phase w ith these changes in
background tem perature. This implies th a t resulting convection cell sizes can be
m anipulated in a controlled m anner by varying the microwave frequency.
We finally tu rn to some observations th a t may be useful for experim ental inves­
tigations.
Given a fluid (such as ethyl alcohol for instance), and the operating
Biot num ber, we now illustrate how Figure 2.7 can be used to predict experim ental
param eters for onset of convection. Consider first a slab w ith fixed thickness d (we
are assum ing here th a t the w idth is much larger than th e thickness). It is useful to
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
27
know th e incident microwave power, Pinc say, required to produce instability. Using
the non-dimensional groups (2.15) and the result Pj„c = q\E$\2 where q is a m aterial
constant, we see th a t R x is proportional to the incident power.
W ith the slab
thickness given, and for known operating microwave frequency, we can determ ine
their ratio, k. Figure 2.7, then, provides a variation between th e required P,nc with
k. Similar conclusions can be drawn from fixing the microwave am plitude |P o | 2 and
the frequency, and using Figure 2.7 to deduce the thickness d required for onset.
Finally we consider possible weakly nonlinear developments and in particular
the emergence of roll or hexagon cellular patterns.
Of relevance is th e work of
Tveitereid and Palm (1976) who consider convection due to an internal heat source
and with boundary conditions of fixed tem perature at one wall the other wall being
insulated. The basic profile is parabolic, then, and the linear stability problem is not
self-adjoint. Due to this, it is shown th a t hexagons are the stable emerging patterns.
It can be shown th a t a sim ilar situation applies here; the analysis is stan d ard and
lengthy and is not included here.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 3
H Y D R O D Y N A M IC S T A B IL IT Y O F A F L U ID L A Y E R :
T E M P E R A T U R E D E P E N D E N T C O M P L E X P E R M IT T IV IT Y
3.1
T he M odel
The physical configuration is th e same as for the constant conductivity case dealt
with in C hapter 2. The fact th a t the dielectric attrib u tes of th e fluid are tem perature
dependent requires scalings of variables th a t are somewhat different from those given
in the constant conductivity case. We consider a viscous incompressible fluid of
constant therm al conductivity kt and viscosity //., bounded between two infinite
parallel plates a distance d apart. The initial tem perature of th e fluid, 60, is assumed
to be the same as the surrounding space, and the corresponding initial fluid density
is denoted by p0. For simplicity, the bounding plates are taken to be transparent to
the impinging microwave field. The geometry is shown schematically in Figure 2.1.
As was argued in C hapter 2 we can average all of th e governing equations
over a microwave period. The resulting equations are a tim e harmonic version of
Maxwell’s equations and th e tim e dependent Navier-Stokes and heat equations. The
la tte r contains the averaged microwave source term (see appendix A).
W ithin this framework, we assume th a t a plane, tim e harm onic electrom agnetic
wave of frequency
and strength E0 impinges normally upon the fluid layer which
fills the region 0 < z < 1. A portion of this wave scatters from the interface
2
= 1,
a portion penetrates the layer and heats th e fluid, and th e remaining portion is
transm itted through the interface z — 0. Let E / denote electric field in the region
above the slab (z > d), E / / denote electric field w ithin the slab
denote electric field below the slab (z <
(0
< z < d) and E ///
0 ).
The electric field generated w ithin the slab acts as a source of energy th a t heats
the fluid. The fluid equations to be solved consist of the Navier-Stokes equations for
heat conducting fluids coupled to a tem perature equation which includes an energy
28
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
29
supply due to the electric field. This electric field is coupled to the heat equation and
can not be solved for exactly. Denote dimensional fluid velocities, tem perature and
pressure by u ' = (u,, u'2, u'3), O' and P ' respectively. Making the usual Boussinesq
approximation, p = /Oo(l —a ( 0 ' —#o)) where a is th e coefficient of therm al expansion
of the fluid, (see Drazin & Reid (1981)) leads to th e following system:
A E '/; - V (V • E '//) + ka [e'(0') + ie"{0')\ E '// = 0,
Du'Po~Dt = ~ 9Po{1 ~ a{9' ~
(3.1)
" V P ' + /ZAU' ’
0 o ))S i3
(3
2)
V • u ' = 0,
(3.3)
DO'
d
d0\
u e "{Q ')^, |2
PoCp D t ~ d x ' j 1dx'j) +
2
|E /; 1 ’
/ 0 .v
(
}
where c is the speed of light in free space, k' = u /c , eo is the perm ittivity of free
space and e! and e" = a / u are the relative perm ittivity and loss factor respectively of
the fluid. If the fluid is purely dipolar, then a is th e im aginary p art of the dielectric
constant. On the other hand, if the fluid is dipolar with conductive losses, then a is
+^
the effective conductivity, Kriegsmann (1993). A =
77.
The boundary
conditions are those of no slip for fluid velocities at the solid walls
u\ = u'2 = u'3 = 0,
z = 0,d,
(3.5)
In the case of a free upper surface th e boundary condition a t z = d becomes
chi',
du'n
du'o
chi',
,
W + a # = 9 ? + #
= 0' “ i = 0'
,
2 = <i’
(3'6)
Convective conditions for the tem perature are imposed a t the walls
flat
k t ^ = h ( e ' - 0 o),
kt^ - = - h ( 0 ' - 0 o),
2
=
0
,
z = d.
(3.7)
(3.8)
In equations (3.2-3.8) above g is the acceleration due to gravity, Cp is th e specific
heat of the fluid and h is the heat transfer coefficient measuring loss of heat due to
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
30
convective cooling. The source term in (3.4) is due to the microwave interaction with
the fluid.
We impose the following continuity conditions on the tangential electric and
magnetic fields a t the slab boundaries; see Ramo et al. (1984), section 3.14.
(V x E '//) x n = (V x E '/) x n,
2
(V x E '//) x n = (V x E '///) x n ,
E '// x n = E '/ x n,
=
1,
2 = 0,
2 = 1,
E '// x n = E '/// x n,
2 = 0,
(3.9)
(3.10)
(3.11)
(3-12)
where n denotes the unit vector normal to th e slab surfaces. The governing equations
are made dimensionless by scaling distances with the slab thickness d, tim e by the
diffusive tim e scale d2/ k where k = kt/(poCp) and the electric field by the incident
am plitude E q . New variables are introduced as
x = x '/d ,
t = Kt'/(d2), u = d u '/«
P = d2/(K2p0)P',
0
= (-l + ^ ),
(3.13)
New electric field variables are introduced as
E = E '/ (e~ik' E q ) ,
k = k'd,
(3.14)
Substitution of (3.13),(3.14) into (3.1-3.12) gives the dimensionless system given
below. The nondimensional groups which enter into the dynamics are the Rayleigh
number, R, the P ran d tl number, P r, th e Biot number, (3, and a nondimensional
measure of the incident microwave power denoted by x- These groups are given by:
R = (agd30o)/(Ki;),
Pt =
u/ k ,
@ = h d / k t,
X = elueo\E0\2d2/ ( 2 k t0o), (3.15)
where eJJ is the dielectric loss factor of th e fluid a t reference tem perature 60. These
scalings give the
following nondimensional system:
A E // - V (V • E //) + k 2 [e'{9) + ie"(0)} E // = 0,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.16)
31
V • u = 0,
(3.17)
r\n
^ - = A 0 + xe" ( 0 ) |E „ |2,
(3.18)
Du
d3
= —V (p + ~^gz) + RPr9\s.+ P r A u ,
Dl = ~ ^ P + n
(V x E //) x n = (V x E /) x n,
2
(V x E //) x n = (V x E ///) x n ,
E // x n = E / x n,
99
Q-Z = W '
— = -(30,
No slip
2
(3.20)
(3.21)
(3.22)
= 0,
(3.23)
Z = 0’
(3’24)
z = 1,
(3.25)
u\ = u2 = u3 = 0,
3.2
= 0,
= 1,
2
E // x n = E /// x n,
2
= 1,
(3.19)
2
= 0 ,1 ,
(3.26)
B asic S ta tes
In th e free space regions 2 > 1 and 2 < 0, the electric field is given by th e real p arts
of
E = E 0 [e -!'(fcz+wr) +
e!(fc2- “r)] j,
z' > d,
(3.27)
E = E 0T e - i(kz+u,T)i
z' < 0,
(3.28)
7
respectively, where E 0 is the strength of the incident field, r is the tim e scale of
the electric field, T is the transm ission coefficient, and 7 is th e reflection coefficient.
B oth T and 7 are to be determined. The electric field which penetrates th e fluid and
interacts with it is given by the real p art of E = [^(z)e~UJT]j, where ip satisfies
iPip
+ k 2 je'(9(z)) + ie"(9(z))j ip = 0,
0<
2
< 1,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.29)
32
Note th a t the averaged microwave source term drives the heat equation. The electric
field from the source is polarized in the y-direction and varies spatially in th e zdirection only; Therefore, the steady-state tem perature 0 depends only on z and the
governing equation 3.1 becomes:
(3.30)
The boundary conditions on ip are derived as in chapter 2 and given as:
(3.31)
(3.32)
dz
Because th e complex perm ittivity is tem perature dependent, the differential equations
governing the electric field and th e tem perature within the slab are now coupled.
As a result, we must solve for the basic states numerically. The base velocity field is
u =
0
. Utilizing the boundary conditions (3.24),(3.25) and (3.31),(3.32), we attain
solutions for ip and 9.
Figures 3.3a-d give plots of th e basic tem perature for different values of k with
unit Biot number and x — -56. Solutions to the above system were obtained by
employing the Newton-Raphson m ethod in conjunction w ith th e ode45 differential
equations solver in M atlab. As can be seen, th e position of the point of maximum
tem perature varies w ith k in the variable conductivity case. D ata points for the
dielectric properties of water, e' and e" , a t various tem peratures were interpolated for
com putational purposes [14]. Figures 3.1 and 3.2 give plots of th e relative dielectric
perm ittivity and the relative dielectric loss factor, respectively, as functions of dimen­
sionless tem perature 6.
In th e case of long microwaves, k «
given in chapter
becomes unity.
2
(section
2 .2 )
1, it can be shown by th e same arguments
th a t the electric field distribution in th e fluid slab
Consequently th e tem perature distribution is governed by the
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Water
76
£6 8
66
60
0.5
1.5
2.5
T e m p era tu re
F ig u re 3.1 Dielectric perm ittivity of water vs. tem perature;
Water
0.5
1.5
2.5
T e m p era tu re
F ig u r e 3.2 Dielectric loss factor of water vs. tem perature;
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
34
Fig b; k = 0.7
Fig a; k = 0.2
0.8
0.8
0.2
0.2
0.1
0.12
0.14
0.16
0.18
dimensionless temperature
Fig c; k = 1.2
0.15
0.2
dimensionless temperature
Fig d; k = 1.7
0.8
0.2
0.06
0.08
dimensionless temperature
0.1
0.05
0.06
0.07
0.08
dimensionless temperature
F ig u r e 3.3 U ndisturbed tem perature profiles for water; /3 = 1.
following equation:
<P6
d? = -* f W '
(3.33)
w ith the boundary conditions (3.24,3.25).
3.2.1
Therm al Runaway
A reasonable question to be asked in speaking about th e stability of a fluid subjected
to microwaves is the following: are th e tem peratures necessary to arouse instability
w ithin the fluid also tem peratures a t which therm al runaway can occur w ithin a
fluid. To arrive a t a reasonable answer, one can investigate the worst case scenario
in which th e fluid is quiescent.
Utilizing equations (3.29) and (3.30) along the
boundary conditions (3.24),(3.25) and (3.31),(3.32) for 6 and ip respectively, we can
produce steady-state curves relating the incident power term x to th e maximum
steady tem perature 0m within the fluid media. Figure 3.4 gives plots of 9m versus
X for k = 1,2,2.5. The fluid considered in this study is w ater since much is known
about its therm odynam ic and dielectric properties. Observe th a t for th e value k = 2,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
35
there exists a power level x beyond which the maximum tem perature jum ps from
the lower branch of each steady-state curve to an upper branch of each curve. The
new value of th e maximum tem perature reached for this curve lies below the boiling
point. For k = 2.5, however, there is a threshold power beyond which the maximum
tem perature will jum p up into the boiling regime. As will be seen in our investi­
gation of the linear stability of the full system (3.16-3.26), these curves indicate th a t
the power levels needed to achieve th e onset of convection fall far below those at
which any therm al runaway phenom ena occurs. If the dielectric perm ittivity e' and
the dielectric loss factor e" are independent of 6 then the curves shown in Figure 3.4
become monotonic. It is the decreasing behavior of e' and e" with 9 which causes
the tem perature distribution to profoundly effect the electric field w ithin th e slab.
This is the root cause of th e non-monotonic behavior of 9m on k. Physically, altering
k can be thought of as altering the thickness of the fluid slab. The neutral curves
were com puted using a secant-shooting method in conjuction w ith a M atlab ordinary
differential equation solver.
An interesting case to investigate is the lim iting case of nearly insulated
boundaries or, equivalently, a small Biot number.
Therm al runaway has been
studied analytically in th is lim iting case in m aterials whose dielectric loss increases
with tem perature (Kriegsmann (1992)).
Studies have also been conducted on
fluids such as w ater where th e dielectric loss and perm ittivity both decrease with
tem perature (Stuerga, More, Lallem ant and Zahreddine (1993)). We reiterate some
analysis here.
In the lim it (/? < < 1), one can construct a leading order asym ptotic approxi­
m ation to th e exact relation between x and 9m. We rewrite x as X ru n fl where /? has
its usual meaning and X ru n = u e " e 0 \ E 0\2d / (2hT0). The system is given:
^ t + k 2 [ e '(6 (z,t)) + ie "(9(z,t))\Tp
=
0,
0
<
*
<
1,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.34)
Steady State Curves
k=1
k=2
k = 2 .5
50
250
200
150
100
300
P o w er (chi)
F ig u r e 3 .4 Steady-state curves for w ater (Bm vs. x); P — 1-
S te a d y S tate C urves
k=1
k=2
k = 2 .5
k=3
100
150
200
250
300
350
Pow er ( c h i j
F ig u r e 3.5 Steady-state curves for w ater ( 6 0 vs. Xrun); P «
1
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Notice th a t the nondimensional heat source of (3.35) contains th e Biot num ber
and hence for ^ ~
0
( 1 ) the source will be very small suggesting th a t the tem perature
will undergo a very slow evolutionary process in tim e t. This suggests the need the
rescale tim e in favor of a slow tim e variable r = (3t. In term s of r , 3.35 becomes:
P ^ = ^
+ X rune"(9)m 2
(3.40)
Formally, 9 and rj) are each expanded in a power series in P [9 = 9q(z , t ) 49i(z, t )/3 -\----- ; xp = tpo + ipiP H
). The coefficients in th e powers of /? are each set
to zero yielding an infinite set of equations. We are interested in investigating the
dynamic behavior of the tem perature to leading order. Substituting th e expansions
above into (3.34),(3.40) and (3.36-3.39) we obtain the leading order system:
^ - + k2 [e'(0o) + it"(Oo)}ip =
0,
0 < z < 1,
(3.41)
(3.42)
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
38
~
2 = 0,
(3.43)
= 0,
2 = 1,
(3.44)
+ ikipo = 0,
dz
dz
0,
=
2 = 0,
— iktpa = —2ik,
(3.45)
z = l,
(3.46)
Notice th a t 0O= ^o(t')- Since 0o does not depend on 2 , the leading order electric field
can be solved exactly. In order to describe the steady-state tem perature to leading
order we must look a t the energy equation and its boundary conditions to th e next
order
0
(/?):
= ^
+
(3.47)
fit)
z = 0,
^ - = 00,
^
= -0 o ,
Integrating (3.47)over th e interval 0 <
2
(3.48)
z = 1,
< 1, we arrive a t ordinary
(3-49)
differential
equation for 0 O or, equivalently, the leading order energy conservation equation:
^
= -20o + w " ( 0 o ) j \ \ E o \ 2d z ,
(3.50)
Since E q is known,this equation along w ith an initial condition on 0o completely
specifies 0o(t). The steady-state equation is:
20o
= Xrun,
e»{00) S l 0\Eo\>dz
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3-51)
39
Figure 3.5 gives steady-state curves for 90 vs. the power term x- T he fluid
studied is water. The curves are given for k = 1,2,2.5,3. Recall th a t the param eter
k can be thought of as the thickness of the slab. Physically we see th a t for k =
1,
the tem perature of the w ater changes continuously with the power term Xrun up to
the boiling tem perature. For k = 2 and k = 2.5 however, the steady-state curve
has an “S” shape. This suggests the existence of a power level beyond which the
leading order w ater tem perature 0 o will jum p from the lower branch of th e neutral
curve to the upper branch of the curve. The tem perature on this new branch is
still below the boiling point for these values of k. For k = 3, however, the upper
branch of the steady curve lies ju st above the boiling point. So, for nearly insulated
boundaries, the model suggests th a t therm al runaway is possible for water. This
is possible because the dielectric attributes create a resonance of the electric field
through their tem perature dependence.
3 .3
L in e a r S ta b ility T h e o r y
The linear stability of th e undisturbed states described above can be analyzed by
adding small perturbations u ', O', P' and E '. Linearizing (3.16-3.26) w ith respect to
prim ed quantities, we obtain:
- A E '„ + V (V • E '„) = - k
^
2
(ec(0)E '„ +
,
+ g t i / = A O' + x ((V’E '/J + ^ * E ',j) e"(0) + O ' ^ f - 1^|2) ,
(3.52)
(3.53)
dtif
*
— = - V p ' - RPrO 'k + P r A u ',
at
(3.54)
V • u ' = 0,
(3.55)
(V x E '„) x n = (V x E ;) x n ,
a = 1,
(3.56)
(V x E //) x n = (V x E'/ 7/ )n,
z = 0,
(3.57)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
40
E 'n x n = E '7 x n,
E /j x n = E 'UI x n,
— =
^
2
2
= -(30',
Rigid —rigid surface
Rigid —free surface
(3.58)
= 0,
(3.59)
= 0,
z =
wf =
= 1,
2
C72
(3.60)
1,
(3.61)
= 0,
- - ■- = 0,
oz1
2
2
= 0,1,
= 1,
(3.62)
(3.63)
where the starred quantities denote complex conjugates and ec = e'(6+9') + e"( 6 + 6 ')
and 9 and ip are the background tem perature and electric field respectively. Taking
the curl of (3.54) twice, we obtain the following equation involving the vertical
velocity component of the fluid and the tem perature:
8
where A i = Jp- +
A w'
= R xP rA x ff + Pr A V ,
dt
(3.64)
This equation along with (3.53) and (3.52)will be used to
determ ine the stability of the system. Note th a t w ith the eigenfunctions w' and
6'
known the continuity equation (3.55) along w ith (3.54) yields the velocity components
u' and v'.
The scalar equations for the electric field components w ithin the slab are:
acv
q2 at
- A E [ - J(z)-g± - G ( z ) ^ = k \(S )E [,
- A E 2' -
~
= k 2 ec(9)E ’2 + k 2 - ( e c(9))9’ip,
- A E ’3 - £ ( J E ' 3) - ^ ( G — ) = k 2 ec(9)E'3,
(3-65)
(3.66)
(3.67)
where J = jg{ec{9))fz {9)/ec-, G = ± ( e c(9))ip/ec
We now consider th e case where E = E (x , 2 , t) in which case u = u (x, 2 , t) and
9 = 9(x, z, t). Note th a t we have dropped the primes from the perturbation terms.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
41
Physically, we can think of convective rolls as corresponding to this two dimensional
construction. From a qualitative standpoint, the two dimensional linear system is
interesting in th a t it describes the m anner in which the perturbation electric field
actually feeds back into the heat equation and affects th e stability of the system
overall. Note th a t for a two dimensional disturbance
and
£3,
£2
is then decoupled from
£1
and is the only component affecting th e energy equation.
9 = T (z ) e iax+at + c.c.,
Given the following ansatz:
£2
=
<£1 (z)etax+<Tt
+
(j)2 {z)e~iax+a~l, w = lV (z)eiax+at + c.c. describing two dimensional disturbances, we
can analyze (3.52),(3.53) and (3.55-3.64) by separation of variables to arrive a t a
linearized ordinary differential system to be solved:
-(D
-(D
2
- a 2 )0! - k 2 ec(9)<l>i = /c2 ^ e c(0)7Y ,
2
- a2W
~ k 2 e * m 2 * = k 2 ±e*(9)Ti>*,
o T + ^ W = (D 2 - a2) T + X(r<t>ie"(9) + ^ V ' ( f l ) ) + T - ^ ( e " ( 9 ) M \
a{D 2 - a2) W = —R P r a 2T + (D 2 - a 2 )2 W,
where D =
(3.68)
(3.69)
(3-70)
(3.71)
The boundary conditions on th e electric field for (k 2 > a2) are:
- i \ J (k 2
- a2 )(j>i = 0,
z = l,
(3.72)
d
<j>1 + i s / ( h 2 - a2 )cj>i = 0,
dz
z = 0,
(3.73)
^ 4 > i+ } / ( * - 0 ) ^ = 0 ,
z = l,
(3.74)
f & - ^ ( a 2 - k 2 )4>x = 0 ,
dz
z = 0,
(3.75)
dz
For (k 2 < a2) they become:
<j)2 * satisfies th e conjugate boundary conditions:
^
2 *+*V(fc2 -
^
2*
aW
= 0,
2
=
1,
(3.76)
- iy/(k 2 - a2 )4>2* = 0,
2
= 0,
(3.77)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
42
for (k 2 > a2) and
^-d>2* + y/(a? =0,
dz
d_
d fa* - y j(a 2 - kz)(f)2 * = 0,
dz
= 1,
(3.78)
z = 0,
(3.79)
2
for (k 2 < a2).
T he boundary conditions for 9 and W are:
W = ^ - W = 0,
dz
6
= 09,
0 = —09,
3.4
2 = 0 , 1,
(3.80)
= 0,
(3.81)
z = l,
(3.82)
2
N um erical Solutions
If the tem perature dependence of the fluid’s dielectric attributes is taken into account
there is very little, analytically, th a t can be done to attain solutions to th e m athe­
m atical problem. In order to solve for th e basic states, we implemented a shooting
m ethod using a Newton-Raphson approach in conjunction w ith th e M atlab ode45
ordinary differential equations solver routine. In the linear stability portion of the
model this approach does work b u t th e system is a tenth order differential system
and the num ber of shooting variables is such th a t the run tim e is too long.
3.4.1
A F inite D ifference Approach
The approach taken here to solve th e linearized problem is to discretize th e system
(3.68-3.82) and pose the problem as a homogeneous system of algebraic equations.
Having constructed the coefficient m atrix, we then im plement N ew ton’s m ethod
to find th e eigenvalue of the problem. For example, the curves of Figure 3.6 are
generated by searching for the roots of th e determ inant of the coefficient m atrix which
is param eterized by the eigenvalue a, all other param eters being fixed. T h e curves
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
43
in Figure 3.7 are generated by shooting on the param eter x-> with a set identically to
zero. Note th a t the Newton iteration procedure is more involved here since x enters
explicitly into the basic state energy equation as well as the coefficient m atrix of the
linearized system.
The following is the discretization for the system (3.68-3.82) on the interval
[0,1] where the discretization step size is A z, Zj — (j —1 )A z and ( j = 1 , . . . , N + 1).
In approxim ating the derivatives, central difference m ethods are employed giving
0 ( A z 2) accuracy on the interior points.
+ 4>ij+i + T jB j — 0,
0 ^ - i + <PhA J +
(3.83)
+ T i B i = °>
(3 -8 4 )
d.9
Tj~i + TjCj -I- Tj+\ + 4>ijDj + <f>2jE j + Wj— — 0,
(3.85)
2 + W j- iF + WjG + Wj+iF + Wj+2 —T jh Aa2R = 0,
(3.86)
where
A j = ( - 2 + A z 2 (k 2 ec(0) - a 2)) , B j = A z 2 k 2 ^ e c{9)il)
Cj = ( - 2 + A z 2(—a - a 2 + ± e » { o m 2 Xj) , Dj = W t f ) ,
F = ( —4 —A z 2 (2a 2 + a / P r ) ) , G =
(6
^
+ A z 2 (4a 2 + 2 a /P r ) + A z 4 (a 4+ a 2 a / P r ))
The boundary conditions of the system are of mixed type and
one-sided difference m ethods give O ( A z ) accuracy.
straight forward
We employ fictitious points
z 0 = —A z and z ^ + 2 = 1 + A z to preserve 0 ( A z 2) accuracy, as explained below.
T he corresponding boundary conditions are for ( j — 1 corresponding to the boundary
z =
0 ):
<j>i 0 =
2 A ziy /k2
— a24>ii + (j>i 2 5
k 2 > = a2,
(3.87)
k 2 > = a2,
(3.88)
0 i o = - 2 A zy/a 2 — /c2 0 Xl + 0 l2,
k 2 < a2,
(3.89)
—2A zy/a 2 — fe2 0 2l + 0 2 2)
k 2 < a2,
(3.90)
(f>20 = —2A z i V k 2 — a 2 0 2 i +
0 2 O=
022’
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
44
wq = W2 , w\ = 0,
(3.91)
To=T2 -2 h p T u
(3.92)
for (j = N + 1 corresponding to the boundary z = 1):
(j>ijv+ 2 = l & z i V k 2 — a2 <j)iN+l + <f>lN ,
k 2 >= a2,
<p2N+2 = - Z & z i 'J k 2 - a2<t>2N+i + ^ 2 jv>
^iw + 2 = ~ 2 A zVa? — k 2 (f)iN+l + <j>iN ,
0 2 jv+ 2
= —2 A z \/a 2 - k 2 (t>*2N + 2 +
# 2 jv>
(3.93)
> = ° 2>
(3.94)
k 2 < a2,
^2
(3.95)
< ° 2>
(3.96)
ww+ 2 = w/v, u;w+i = 0,
(3.97)
T n + 2 = Tn — 2h(3Tw+i,
(3.98)
The unknowns 0 io, <£2 0 >wt>’ ^ 0 and
0
iat+2 > ^ 2 n+ 2 ->wn+ 2 , T n +2 can all be eliminated
from (3.87-3.98) by employing central differences of th e equations (3.83-3.86) at z = 0
and z = 1 respectively. The modified equations are as follows:
for j =
1
the system is given as:
+ 20i 2 + T \B i = 0,
(3.99)
4>'2 l A ? + 2(j>*22 + T iB l = 0,
(3.100)
TxC; + 2 T2 + 0 u D i + ^ E j . = 0,
(3.101)
Since w\ = 0, th e momentum equation does not enter here.
For j = 2 the only
equation modified in th e system (3.83-3.86) is (3.86) which becomes:
W2 ^G + 1) + W3 F + 1U4 —T^Az^o^/? = 0,
(3.102)
for j = N equation (3.86) becomes:
ujyv—2 4" Wff—i F + Wff(G + 1) —Tff Az^o^TZ = 0,
(3.103)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
45
F igure 3.6 Growth rate curves, rigid-rigid case; /? = 1; k = 1; P r = 7.
for j — N + 1 the system becomes:
2 0 1
2 2 *
+
N
+
^ lw + i-^ w + i
+
T a t+ iB a t+ i
=
2 0 2 jv +
0 2 w + i - ^ jv + i
+
? V + i-5 w + i
=
T
n
+ i C 'n + 1 + 4>i n + 1D
n
+i +
0,
0 2 //+ i-E W
( 3 .1 0 4 )
( 3 .1 0 5 )
i
=
0 ,
( 3 .1 0 6 )
where A' = Aj + 2 i A z \ / k 2 — a 2 if k 2 > a2, A'- = A j — 2A z \ / a 2 — k 2 if k 2 < a2,
Cj = Cj —2 A z/3. The discretized system is now 0 ( A z 2) accurate.
3.5
R esu lts
The stability characteristics were obtained by fixing the physical param eters
including the microwave power which is represented as the nondimensional variable
X ■ For different wavenumbers a, the eigenvalue a, which is complex in general,
was com puted by iteration. O ur numerical results produce real values of a alone,
indicating th a t there is no propagation of disturbances down the slab. This would be
expected even with th e introduction of tem perature dependent dielectric attributes
since there is no physical mechanism present to motive propagation of a disturbance
in a preferred horizontal direction. As a result, neutral curves can be traced by
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
46
specifying a = 0. A typical set of results is given in Figure 3.6 for different power
settings (values of x) and physical param eters corresponding to water. W ater was
chosen to be th e subject fluid as there exist d a ta for the dielecric p erm ittivity and
loss factors (e1 and e") for 0 < 9 < 4 for a microwave frequency of around 2.45
GHZ. As can be seen, larger values of x support a larger band of unstable modes
which narrows as the value of x decreases-see th e curves corresponding to x = -58
and x = -55. As the microwave power is decreased further, it reaches a value below
which no mode becomes excited (slightly below x = -5).
Next we use typical stability plots, as in Figure 3.6, to present neutral curves
which have a = 0. For a fixed x if is seen th a t there are a t m ost two values of
the wavenumber where the flow is neutrally stable, and below a critical value of
X,
denoted by X o the flow becomes linearly stable to all wave numbers. N eutral
curves of x against a for different Biot numbers (other param eters being fixed) are
considered first. In Figure 3.7 we show neutral curves for /? = 0.1 and 1.0 when the
non-dimensional param eter k = 1.0. The unstable and stable regions are located
above and below the given neutral curve respectively. In fact each neutral curve
provides a critical minimum point (ac,x c); flows corresponding to Rayleigh numbers
below Xc are stable. T he variation of (ac, Xc) in param eter studies provides a lot of
inform ation about the stability characteristics and is considered later. Physically, /?
may be associated w ith the rate at which heat is convected away from th e slab and
hence an increase in the Biot num ber has a stabilizing effect on the system . T h a t is,
for a given value of x th e window of unstable modes contributing to th e disturbance
in the fluid decreases.
Figure 3.8 gives th e variation of the critical values of Xc w ith k, th e ratio of
slab thickness to microwave length. Figure 3.9 is a plot of the critical wavenumber ac
versus k, and Figure 3.10 presents the variation of the height where the background
tem perature achieves its maximum w ith k. The critical values for x are those values
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
47
Biot num beral
Biot nu m b ers.1
F ig u r e 3 .7 Neutral stability curves in the x~a plane for water; k = 1.
at which the onset of instability is achieved and exactly one perturbation mode
corresponding to the critical wave number ac becomes excited. We note th a t the
plots in Figure 3.8 have an oscillatory like behavior as k varies. The plots where
constructed based on physical param eters for pure water. We com pare the plots
obtained when water is assumed to have a constant dielectric conductivity w ith those
obtained when tem perature dependence is taken into account. Note th a t in the case
of constant perm ittivity the values used are at the reference tem perature 90 = 0.
(see Figures 3.2 and 3.1). Observe th a t for certain values of k
( .2
< k < .27 and
.43 < k < .5) the assumption of a constant complex dielectric perm ittivity leads to
lower estim ates on of the power term x a t which the onset of instability occurs th an
does the assum ption of a tem perature dependent complex perm ittivity. For all other
values of k, the opposite is true. The two curves are out of phase for sm aller values
of k, b u t become more syncronized as k gets larger. This may be explained by the
fact th a t th e largest fluctuations in the maximum tem perature Tm a t onset are seen
for smaller values of k (see Figure 3.11) and hence, the dielectric p erm ittivity e'(9)
and the dielectric loss factor (."(9) vary significantly from their values a t th e reference
tem perature. As k increases, however, these fluctuations decrease as does the overall
maximum tem perature of the fluid as seen in the figure. As was mentioned in section
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
48
variable dielectric permittivity
constant dielectric permittivity
F ig u re 3.8 N eutral curves in the Xc ~ k plane for water, rigid-rigid case: f3 = 1;
k = 1.
3.2.1, there is a strong coupling between the tem perature 9 and the electric field
through the fluids dielectric properties e'(9) and e"(9). These dielectric properties
have a significant effect on the resonance of ip within the fluid layer and hence, on
the shape of the basic tem perature profile. For example, for .3 < k < .4, it is seen
in Figure 3.10 th a t when the tem perature dependence of the dielectric properties of
water is taken into account, the height of maximum tem perature within the fluid layer
is less than if this tem perature dependence is neglected. If the height of the maximum
tem perature is lower, then a larger unstable fluid layer resides at th e upper boundary.
As a result, less power is necessary to achieve an onset of instability (see Figure 3.8)
and such an instability can be achieved at lower tem peratures. The fact th a t the
curves of Figures 3.8-3.11 are in phase for both the case of constant complex perm it­
tivity and tem perature dependent complex perm ittivity, is justified by arguments
similar the those given for the constant conductivity case (see chapter 2 ).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
49
variable dielectric permittiyit)
constant dielectric permittivn
*
I
F ig u r e 3.9 Neutral curves for w ater (critical wave number vs. k), rigid-rigid case:
P
=
1 ; jfe = 1 .
variable dielectric permittivity
c o nstant dielectric permrttivit
F ig u r e 3.10 Neutral curves for w ater (height of maximum tem perature vs.
rigid-rigid case: /? = 1 ; k = 1 .
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
k),
50
variable dielectric permittivit
constant dielectric permittivi
2
o.<
F ig u r e 3.11 N eutral curves for w ater (maximum tem perature vs. k), rigid-rigid
case: (3 = 1; k = 1 .
3.6
D isc u ssio n
The model considered in this chapter takes into account more of the physics involved
in the interaction of a fluid w ith microwaves in th a t it incorporates th e tem perature
dependence of the dielectric attrib u tes of the fluid. Such a model is of im portance
in order to achieve greater accuracy on the power requirements to achieve the onset
of instability where tem perature variations within the slab are significant. As was
seen for water, there are intervals of k where the monotonic decreasing behavior of
the dielectric loss factor e"(Q) and the dielectric perm ittivity e'(9) w ith tem perature
reduces th e critical power level Xc necessary to achieve an onset of instability. So,
in these regimes, ignoring the tem perature dependence of the dielectric properties of
w ater on tem perature would lead one to over estim ate th e incident critical microwave
power necessary to achieve onset. There are, yet, other intervals where assum ption of
constant dielectric properties would lead us to underestim ate this critical power. In
carrying out industrial applications where one often wishes to use th e least am ount
of power possible, the results discussed have im portant implications.
Figures 3.4 and 3.8 indicate th a t the power levels needed to achieve an onset of
convection fall far below those necessary to experience therm al runaway phenom ena
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
51
in water. W hat is interesting about Figures 3.4 and 3.5 is th a t they illustrate the
im portance of geom etry in achieving therm al runaway. I t was seen th a t by varying
the thickness k of th e fluid layer, one can alter bo th th e tem perature and microwave
power level at which therm al runaway takes place. This m ay have im portant im pli­
cations in commercial processes such as microwave processing of foods or other waterbased materials.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 4
A M ODEL FO R M ELTING OF SOLIDS U SIN G M ICROW AVES
4.1
A P hysical O verview (T he Stefan C ondition)
The physical system presented here consists of a two dimensional m aterial of infinite
extent in the x direction which is heated by a microwave source. Specifically we
assume the existence of a tim e harm onic plane wave which impinges normally to the
top face of the m aterial. This electric field is taken to be polarized in th e x direction.
The m aterial consists of two phases, liquid and solid, of a substance. The interface
bounding the two phases is referred to as the m elting front. One m ain objective
is to track the position of the melting front in tim e. As can be seen in figure 4.1,
there exist energy equations for both the solid and liquid phases which govern the
tem perature distribution in each of the respective regions. We let <j>i and <f>2 represent
the electric field distributions in the liquid and solid phases respectively, and T\ (z)
and T 2 ( z ) are the corresponding tem perature distributions in the respective regions
and Ki and k 2 the corresponding therm al conductivities. Figure 4.2 illustrates the
physics governing the movement of the m elting front S(t). At the melting front there
exists an energy flux into the front from the liquid side —
out of the front on the solid side
and an energy flux
If these fluxes are not in balance then the
excess energy is consumed in facilitating a phase change o f th e solid into liquid a t a
rate p
a
Hence the difference in these fluxes a t tim e t determines th e velocity of
the melting front
This is th e Stefan condition:
dTi
dT2
dS
- K i — + K2— = p a — ,
dz
dz
dt
, .
(4.1)
a is the latent heat o f m elting and p is the density of th e solid. Note th a t the
front m aintains a constant tem perature known as th e freezing tem perature which is
denoted as Tj.
52
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Z=Or
T1a.z).«la.z)
Melted region (I)
t 2(u).$2(u)
Frozen Region (II)
z=S(t)_
z=d
T j(t,z ) . T ,(i,z )-T em [> e n itu re d istrib u tio n s
<J>j(t.z) <]>2 ( ttz )-E ]c c tric field s
F ig u re 4.1 The geometry of the melting problem.
p a d s-e n e rg y consumed to melt solid
Kf*"l z -energy flux into region
K2 T2 z -energy flux out o f region
- k ;T1 z +
p a S j — The latent heat balance (Stefan) condition
p - density o f solid
a - latent heat o f melting
F ig u re 4.2 The physics governing the moving front.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
54
4.2
T he M athem atical M odel
We present the dimensional model governing the electric fields and tem perature
profiles w ithin the liquid and solid phases of th e material.
+ k 2 (a i e\(T{) + c n t'K T D M = 0,
dz2
dP(j}'2
+ k 2 (a 2 e'2 (T ') + a 2 4 ( T ^ ' 2 = 0,
dz 2
=0
+^ w
0
<
2
< S '(i),
S '( t ) <
2
< 1,
2’ °<*<*'«>•
+
« • ( * ) < * < 1.
(4.2)
(4.3)
<«>
(«)
where e\ and e'{ are the relative dielectric perm ittivity and loss factor, respectively, of
the liquid region and e2' and e2 are the relative dielectric perm ittivity and loss factor,
respectively, of the solid region: Also, a i and a 2 are values of the dielectric p erm it­
tivities of these respective regions a t th e reference tem perature which wedenote as
T0,
and o\ and a 2 are the values of th e loss factors a t this tem perature; k = ui/c
where to is the frequency of the im pinging electromagnetic wave and c is the speed
of light in free space; ci and c2 are the specific heats of the liquid and solid regions
respectively and e0 is the dielectric p erm ittivity of free space. Note th a t (4.4),(4.5)
are driven by tim e averaged electric field sources. The tim e averaging is justified b)r
the fact th a t the plane wave oscillates on a time scale which is much smaller th a n the
therm al diffusive tim e scale of both th e solid and liquid phases of th e dielectric (see
Appendix A). M anipulating Maxwell’s equations, we arrive a t ordinary differential
equations governing the electric fields in the solid and liquid regions.
The electric field in the free space above th e m aterial has th e form:
(j>0 = E 0 (eikz + 'ye~ikz)
(4.6)
The transm itted electric field in the free space region below the m aterial has the
form:
4>o = E 0 T eikz,
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
(4.7)
55
where E 0 is the am plitude of th e incident wave,
7
is the reflection coefficient and
T is the transm ission coefficient. A t dielectric interfaces, we impose conditions of
continuity of both th e tangential electric and m agnetic fields (see Van Duzer, Ramo,
section 3.14). C ontinuity of the tangential m agnetic field is equivalent to continuity
of th e norm al derivative of the electric field a t z = 0 and z = d. These conditions
yield the following equations:
(j>\ = E 0 (eikz + j e ~ ikz)
dz
Elim inating
7
= ik E 0 (eikz -
7
z = 0,
z = 0,
e - ifcz)
(4.8)
(4.9)
in favor of <j>[, the boundary condition 4.9 becomes:
^ + ik<t>\ = 2ik E 0,
dz
z = 0,
(4.10)
Applying these sim ilar conditions at the bottom face gives:
^ - i k < j > ' = 0,
z = d,
(4.11)
On the interface S '( t) separating the liquid and solid phases, the boundary conditions
are those of continuity of electric and magnetic fields, ie
« > ',= & ,
^
=
z = S'(*).
(4-12)
The boundary conditions on th e tem peratures T[ and T \ in the liquid and solid
portions, respectively, of the m aterial are given by:
FfT'
Kl - ^ = - h { T [ - T 0),
T i '= T ' = 0,
T ' = Tbw,
dT[
dT^
- « ! — + k2— =
z = 0,
(4.13)
z = S'(t),
(4.14)
z = d,
(4.15)
dS'
,
z = S {t),
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4.16)
56
where h is the heat transfer coefficient a t the upper surface and TbW is some fixed
tem perature.
The following scalings are introduced:
z = z '/ d
t = t'K2/(<Pp2c2)
Ti = (T{ - T f ) / ( T 0 - Tf)
fa = f t / i E o )
fa = f a / ( E 0)
T 2 = (V2 - T f ) /( T 0 - Tf )
S = o2/(Ji
A = k 2/ k \
p = PiCi/(p2c2)
P = hd/Ki
L = (a \ ) / ( ( T 0 - Tf )c2)
f = a 2/ a i
P = ai<jje0 cP\E0 \2 / ( 2 (T0 - Ts )tcx)
W ith these scalings, we have the following dimensionless system:
+ ia ^ T ^ fa =
^
0
,
0
d?fa
+ k 2 (a\€e' 2 {T2) + ic i 8 e2 {pT2))fa = 0,
dz 2
<
2
< S(t),
S'(t) < z < 1,
(4.17)
(4.18)
= ^ + P e ' l ( T 1 )\fa\2,
0 < z< S (t),
(4.19)
dT 2
d 2T2
5
.. |2
-gj- = -q ~£ + P j h i f a n f a ] ,
. .
S (t) < z < 1,
(4.20)
p \ ^
z = 0,
(4.21)
z = S(t),
(4.22)
z = 1>
(4.23)
~ r = P(Tr - 1),
fa — T 2 = 0,
T2 = Tbw,
dz
+ ik fa = 2 ik,
fa = fa,
dfa
0,
. .
3 “ s (t)0
,
2
=
(4.24)
(4.25)
z - S(t),
dfa
^ P - ik fa =
dz
z=
1,
(4.26)
(4.27)
The nondimensional Stefan condition is:
dT\
dfa
dS
-~ ^ +^
= LT t’
z = 5 ( t) ’
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4.28)
57
A t tim e t = 0, we let Ti = f ( z ) and T 2 = g(z). This system comprises th e full initial
value problem to be solved.
4.3
4.3.1
T he C ase o f C onstant D ielectric P e r m ittiv ity
S tead y-State Solutions
It is of im portance in m elting applications, to know how much of a m aterial will m elt
for a given power level of a microwave source. If there exists only one m elting front
within a
m aterial, it isenough to determ ine the steady-state position of the m elting
front Seq fora given power level P. To find the steady-state solutions of the system
(4.17-4.28), we set the tim e derivatives of this system to zero. The equations are
given as:
—j - ^ - + k 2 (aie '1 + icrl e")(l>ieq = 0,
^ z 2~ 4" k 2(a l€ e2 +
d^2)<l>2eq =
0 < z < S eq,
0)
&eq < Z < 1,
< S ,„
(4.29)
(4.30)
= - P « f I#, « / .
0 <
= - P jt'ilfo ,/
S „ < Z < 1.
(4.32)
z = 0,
(4.33)
2
(4.31)
p fp
^
l
= P(Tleq- 1),
T \eq = I 2 eq = 0,
z = S eq,
^ 2 eg = Tbw,
d<t>leg
+ ikfii
= 2ik,
<l>lea
1eq = <f>2ea
Y2eq
dz
z = 1,
z = 0,
(4-34)
(4.35)
(4.36)
«* = S*■Jeqi
eq,
(4.37)
*= s
(4.38)
dz
- ik<j)2eq = 0,
z = l,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4.39)
58
The nondimensional Stefan condition is:
-
^
+ A ^ 2 = °,
z = S„,
W ith th e assum ption of constant dielectric attributes,
(4.40)
e2, e",and e2 are unity.
Observe th a t equations (4.29) and (4.30) are decoupled from the heat equations
in the solid and liquid regions. Using these, together w ith the boundary conditions
(4.36)-(4.39), we can deduce an exact solution for the electric fields in both the liquid
and solid portions of the material. The solutions are given as:
<t>Uq = c 1e - ir' z + c2 eir' z,
(4.41)
</>2eq = c3 e - ir'* + c4 eir**,
(4.42)
where ci through c4 satisfy the following linear algebraic system which is derived
from equations (4.37-4.39):
( Ti + k
—r i + k
0
gir i Seq
g—iT1Seq
_gir'2 Seq
r v r >s-
- r 1 e - ir‘s««
- r 2 eir*s'<
o
0
( - k + r 2 )e ir 2
V
fc A
/ 2k\
_Seq
C2
0
r 2 e -iV*s"
c3
0
0
1
- { k + r 2 )e~ir>) \ ^4 /
lo y
where Tj = y/ai + i a i and r 2 = V ^ai +
T Uq(z) = - Pe'[F{z) + A z + B ,
(4.43)
T 2 eq(z) = - P j t ' l G ( z ) + C z + D,
(4.44)
where
|c i| 2 e-2/m(r i )z
Re{c\c2) cos(2i?e(Fx)z)
4 / m ( r 1) 2
2Re{Y1y
Im {c\c2) sin( 2 i 2 e ( r 1 )z)
|c2 |2 e2/7n(r i >*
2 i 2e ( r j ) 2
+ 4 J m ( r 1) 2
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
59
r( n
|c3|2e 2/m(r =)2
-Re(cgc(4)) cos(2i?e(r2)z)
2JR e(r2)2
Im{c\c{A)) sin(2i?e(r2)z)
|c4 |2e2/m^r2^
4/m(r2)2
2i?e(r2)2
,4
'4 ~
—0 ( P F (0)+ 1)+ /?P F (S cq )+-P ^r (0)
0 S e, + l
^
_
P<5/A(G(1)—G(o))+7),,i,
°
~
I —Sen
4/m(r2)2
„
B ~
n
^
Starred quantities denote complex conjugates.
P F ( 5 e<, ) - S e, P ^ ( 0 ) + S e,/? (P F (0 )+ l)
0 S eg+ 1
Ptf/A (-C 7(l)5<;„ + G (5 eo) ) - r tu,Seo
1- S
~
e n
W ith these solutions known, the
steady-state Stefan condition (4.40) becomes
~ { P ^ + A) + 5 P ^ + \ C = 0,
dz
dz
(4.45)
Equation (4.45) constitutes an eigenvalue relation relating the steady front position
S eq to th e power term P . We can construct steady-state curves relating the steady
front position S eq to the microwave power term P using Newton’s method. Specif­
ically, we iterate through values of S eq ranging from .02 to .9 and, for each iteration,
perform a Newton’s iteration on the power level P . Figure 4.3 gives the steadysta te positions of the melting front denoted by S eq as a function of power produced
by the microwave source for three different Biot numbers. Recall th a t T/,w = —.5
represents the value of th e tem perature T2 at z = 1. The param eter values used in
our numerical simulations are those corresponding to w ater and ice: a \ = 80,
6
=
.01, £ = .04, A = 4, fi = 2. The microwave frequency of consideration is about 2.45
GHZ. It is interesting to note the oscillatory nature of the curves, particularly the
curve corresponding to Biot number 0 = 1. This curve indicates th a t for certain
power levels, there exist multiple steady-state positions for the m elting front. As the
Biot num ber is increased, the steady-state curves become more monotonic. W hat is
interesting here is th a t these curves all seem to intersect near a particular value of
S eq. In fact, this value SpiVOt can be calculated by considering th e case where P = 0,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
60
|
3000
F ig u re 4.3 Steady-state curves (constant complex perm ittivity with varying Biot
number): Ttw = —.5; k = 1.
which can be thought of as the absence of a microwave source, and by letting j3 go to
infinity or, equivalently, applying a Dirichlet condition on T\ a t z = 0 (ie. T\ = 1).
The steady-state position of the front for these param eter values is exactly SPiVOt.
Notice th a t for steady-state front positions S eq < Spivot the power level needed to
achieve a given steady-state position decreases w ith increasing Biot number whereas
the power level increases for steady-state positions S eq > SPiVOt . This is due to the fact
th a t at z =
0,
th e boundary condition can assist in the propagation of the melting
front up to the value S eq = SPiVOt so an increase in the Biot num ber /? enhances th e
propagation. For values S eq > SPiVOt, an increase in the value of the Biot num ber f3
can only inhibit the m elting front propagation.
Figure 4.4 gives steady curves of S eq vs. P for different boundary values
of the tem perature Ti- As the boundary tem perature decreases, the power levels
needed to achieve a steady-state value S eq increase. This is to be expected as th e
heat flux across the m elting front is increased which inhibits th e movement of this
front. The steady-state curves become less monotonic as th e tem perature a t th e
boundary is decreased.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
61
{
600
clMdyfrontportion(S^
F ig u r e 4 .4 Steady-state curves (constant complex p erm ittivity w ith varying
boundary condition): (3 = 1; k = 1 .
dielectric interface (z= i)
perfectly conducting interface (z*1)
i
■toadyfrontposition(S^)
F ig u r e 4 .5 Steady-state curves (constant complex perm ittivity cases): f3 = 1; k — 1;
T bw = - .5 .
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
62
Figure 4.5 compares the steady-state curve for a m aterial layer whose bottom
interface is a dielectric interface w ith a m aterial whose bottom interface is perfectly
conducting.
surface.
Physically, one might think of a m aterial layer resting on a m etal
M athematically, the presence of a perfectly conducting surface
can be modeled by changing the boundary condition (4.27) to <j>2 — 0 a t
2
(2
=
1)
=
1.
As
can be seen in Figure 4.5, the introduction of a perfectly conducting surface has a
noticeable effect on the power necessary to achieve a desired front position S eq.
4.3.2
A sym p totic Lim it o f Large Stefan N um ber
In the case of a Large Stefan number L, we can introduce a small param eter (e =
1/
L)
into the system (4.17-4.28). As will be seen, approxim ate solutions can be found
analytically for the electric fields and the tem perature distributions of both the liquid
and solid regions of the slab. Consequently, the Stefan condition, at leading order,
translates into an ordinary differential equation governing th e front position which
can be solved using standard numerical techniques. The assum ption th a t L is large
m otivates th e scaling r = et. S ubstituting this scaling into the system (4.17-4.28)
gives the following system:
0<
2
< S(t)
S( t) < 2 < 1,
(4.46)
(4.47)
(4.48)
(4.49)
(4.50)
Ti = T 2 = 0,
2
= S(t),
(4.51)
T 2 = Tbw,
2
—1,
(4.52)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
63
dz
+ ikcfii = 2ik,
^i =
d<t>\
02,
2
2
= 0,
(4.53)
= S (r),
d(j)2
(4.54)
_, .
dT= 17’
,.
z = S(T)’
(4'55)
z = 1,
(4.56)
^
- i k * = 0,
dz
dTi
dT 2
dS
, *
=
* = 5 ( t ) -
,
.
( 4 '5 7 )
The problem is now posed as one with a small param eter e. We now solve the system
by regular perturbation theory. We expand the variables of the system in powers of
e as follows:
T\{ z , t ) = T ° ( z , t ) + e T i ( z , r ) + 0 ( e 2),
( 4 .5 8 )
T2(z , t ) = T2°( 2 , r ) + e T 21 ( 2 , r ) + 0 ( e 2),
( 4 .5 9 )
<f>i{z,r) =(f>°l (z ,r) + e(j)\{z,T) +
0 ( e 2),
( 4 .6 0 )
+etf>l(z,T) +
0 ( e 2),
( 4 .6 1 )
<£2 ( 2 , 7- )
= # ) ( 2 , t )
5 ( r ) = S 0(t ) + eSi(r) + 0 ( e 2),
( 4 .6 2 )
S ubstituting (4.58-4.62) into system (4.46-4.57) and equating like powers of e we
a ttain th e successive systems to be solved. Note th a t functions are represented by
a Taylor series in
2
about S 0 a t the boundary
2
= S . In general, a function f ( z ) is
expressed as f ( z ) = f ( S 0) + %{S0){z - S 0) + j £ { S 0){z - S 0 ) 2 + 0 ( z - S 0 ) 3 a t z=S.
This enables us to deduce boundary conditions at each order of e. To leading order
the system is as follows:
dz2
+ k2{ai + icri)4>\ = 0 ,
<P<t>l + k 2 (ai£, + ia\5)(j>2 = 0,
dz 2
d 2 T° _
= -P|<£?|2, 0 <
dz2
0<
2
2
< S o(t ),
(4.63)
S 0 (r) < z < 1,
(4.64)
<
(4.65)
5 o(t ),
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
64
d ”? 1®
dz2
r ,^
JirrpQst
,0 |2
S q(t ) < z < 1,
= - P - e " ( r 2° ) |^ |2,
A
dT°
= 0 (7 ° - 1),
dz
T ° = r 2° = 0,
z = 0,
(4.67)
z = So(r),
(4.68)
z = 1,
(4.69)
T ° = T - bw,
^
+ ik(j>\ = 2ik,
dz
4>\ = 4>°2,
z = 0,
(4.70)
z = S q{t ),
d(j>2
^ 1
(4.66)
(4.71)
___ „
2 -
=
/ , -n\
, v
s»W ’
^
- i k $ = 0,
dz
dTf
dT°
dS 0
1 4- A—^
,
dz
dz
dr
2=
<4 - ' 2 >
(4.73)
1,
_
z = S0,
(4.^4)
The 0 (e) system is as follows:
dz 2
dz 2
+ k 2 (a\ + ioi)4>\ = 0,
0 < z < 5o(r),
(4.75)
+ k 2 (ai£ 4- ia\5)(l>l = 0,
■S’o (r) < z <
(4.76)
= l 3 ~+ 2PR e^ ? ^ * ) ’
^
= ^
+ 2 P ^ R e (^ * ),
= pTl
^
r1T°
1,
0 < z < 5 °(r )>
(4 ‘77)
5 0 (r) < z < 1,
(4.78)
z = 0,
(4.79)
dT°
T l ( S 0) + ^ ( 5 o ) 5 i = ^ ( S o ) + - ^ ■ ( S 0 )S 1 =
T2 = 0 ,
az
4-
0,
z — 1,
= 0,
(4.81)
z = 0,
(So) 4- ^ ■ ( S 0 )S 1 = <f>l(So) + ^ { S
(4.80)
(4.82)
0 )S U
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4.83)
65
M (So) + 2^
dz
(So)Sl = M ( < y +
—ik(j)\ = 0,
2
2!M
( S 0)S l,
= 1,
(4.84)
(4.85)
. \dT«
.
d?T2 rr, 1
d5i
dTl1 ( 5 0 ) 4 - 2 5 x ^ 5 ( 5 0 ) + A ——(So) + 25i , „ (So) =
dz
dz 2
.t
(4.86)
S tarred quantities denote complex conjugates and Re denotes the real p art of the
quantity on which it operates. We now restrict our attention to the leading order
system, keeping in mind th a t th e higher order systems can be solved to achieve higher
order accuracy of the solutions.
Observe th a t the electric fields <j>\ and 4>%along w ith tem perature distributions
T ° and T2°, of the leading order system, satisfy the same equations as the steadystate system (4.29-4.39). Hence, <j>\ = <j>Uq, <f% = 4>2eq, T? = T leq, and T2° = T2eq.
<$>) T ° and T ° depend
The system (4.63-4.74) is a quasi-static system in th a t
param etrically on the front position S 0( t ) . Since the tem perature and electric fields
of both regions are known functions of S q ( t ) , th e Stefan condition (4.74) becomes the
governing ordinary differential equation for the leading order melting front trajectory
- ( P ^ + A) + 5 P ^ - + XC = ^ ,
dz
dz
dr
(4.87)
E quation (4.87) and an initial condition on S 0 constitute an initial value problem to
be solved. Solutions are found using a M atlab ordinary differential equations solver.
Figure 4.6 gives plots for the position of the melting front S{t) vs. its rate of
change S'(t) in th e lim it of large Stefan num ber (L »
1 ).
The plots are given for
two different power levels. As can be seen for P = 400, there exist three steady-state
positions for the melting front. The equilibrium points S = .15 and S = .47 are
stable nodes while S = .21 is an unstable node. If the power P is increased by 40
units the num ber of equilibrium points is reduced to one. This point S = .48 is a
stable node. So, for a relatively small change in the power P , the dynam ical nature
of the system (4.104-4.115) is altered significantly.
Figure 4.10 shows the front trajectory in tim e for L »
1.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
66
front S{t)
Figure 4.6 M elting front position S vs. derivative S': (3 = 1; k =
4.4
4.4.1
1.
T he Case o f Tem perature D ependent D ielectric A ttrib u tes
S tead y-State Solutions
The calculation of the steady-state solutions in this case where the dielectric
attributes of a m aterial depend on tem perature, is more com plicated th an the
calculation for the case of constant dielectric properties. This is due to th e fact
th a t th e tem perature fields and th e electric fields of the liquid and solid regions
are now coupled together. We will use the same notation to denote th e equilibrium
solutions as we did in the previous section. Taking the tim e derivatives of th e system
(4.17-4.28) to be
zero, we arrive a t the steady-state system:
+ * 2 ( a xe i ( r le,) + « 7 lC?(Tle,) ) * le, = 0,
^ 3
+ k?(a 1t f 2 (T2eq) + za 1 (5e"(T2 e9 ))0 2e7 = 0,
rPT
- ^ 3
rPT
if
= - P ^ ( T 2eq)\cj>2eq\ \
Tleq ~ Tzeq = 0)
2
< S eq,
S eq < z <
1,
0 < z < S eq,
= - P e " ( T le9 )|0 le/ ,
^ ~ 3 _ p (Tl _ i),
0<
2
(4.88)
(4.89)
(4.90)
S eq < z < 1,
(4.91)
= 0,
(4.92)
z =
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4.93)
67
constant conductivity
variable conductivity
* tsoo
flMdy(rentposition(S^)
F ig u re 4 .7 Steady-state curves (constant and variable complex perm ittivity cases):
(3 = 1; k = 1; Tbw = -.5 .
2
= 1,
2
+ ik 4>\eq = 2 ik,
(filea
^ 2 eo)
(4.94)
= 0,
^ = s .eqi
d(pleq _ d fc eq
dz
dz
d&2 eq
- ik(j>2 ea = o,
dz
2
— S eq,
(4.95)
(4.96)
(4.97)
2
= 1,
(4.98)
2
— §eqi
(4.99)
The nondimensional Stefan condition is:
d ju q + ^ dT2eq
—0 ,
dz
dz
The problem posed here is an eigenvalue problem relating the steady-state front
position Seq to the power term P . T he steady-state tem peratures, electric fields
and melting front position are found by implementing a shooting m ethod on the
system (4.88-4.99). Specifically guesses are made for (f>ieq(0), 7 \e9 (0) and the power
P . So, for a given value of the param eter S eq the shooting m ethod converges to a
certain value of P . To solve this system we employ a Newton Raphson m ethod in
conjunction with a M atlab ODE solver routine.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
68
Figure 4.7 displays the steady-state curve with TbW — —.5 for the cases of
tem perature dependent dielectric perm ittivity and constant dielectric perm ittivity.
The nonmonotonic behavior of both curves is very pronounced. There is a significant
difference in th e steady-state power levels predicted by the curves, which demon­
strates th e im portance of including the tem perature dependence of th e complex
dielectric perm ittivity in the m athem atical model.
4.5
A Fixed Front M ethod For Tracking th e M oving B oundary
The main idea of the front fixing m ethod is to make a transform ation of variables
in order to transform the moving boundary value problem (4.17-4.28) into a fixed
boundary value problem where th e boundary S(t) will now enter into the partial
differential equations of the system. For S ( t ) < z < 1, we let £ = (z — S(t ) ) / ( 1 —
S(t))). For 0 < z < S(t), we let rj = z / S( t ) . W ith these transform ations, th e region
S(t) < z < 1 is m apped into 0 < £ < 1, and 0 < z < S(t) is m apped into 0 < ij < 1 .
The melting front is a t C = 0,
derivatives in t and
2
77
= 1. The transform ation of the tim e and space
are as follows:
d
d
dS Q—I d
.
at~*Ft+ d t T ^ s d ( ’
.
...
0<c<1'
(4100)
0<f<1>
l - s - i i s ?
£ -> 3 3 ?
.
(4101)
0 < ” <1- .
(4102)
0<4<1’
(4103)
W ith these transform ations the system (4.17-4.28) becomes:
drj2
4_ p 2 1.2 /
+
+ u r^ T i)^ ! =0,
d?<j>2
+ (1 - S ) 2 k 2 {a^e' 2 {T2) +
dQ2
= 0,
0<
< 1,
(4.104)
0 < C < 1,
(4.105)
77
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
69
,&rx d s \ j r r Xs_
~
dT 2
dt
l a 2^ , n , „ ^ , IJL|2
+
=
dSC -ldT2 _
d t l - S dC
(1
^
1
d 2T 2
- S ) 2 dQ2
0
ps
A (
= 5/?(Tx - 1),
=
,
0
, r
< c <
(4-106>
, 1
1,
0,
(4.109)
C = 1,
^
+ i k S f a M = 2ikS,
dr)
4>i(*i) = M Q ,
C=
(1 -5 )^ = 5 ^ ,
>2
=
0
rfTi 1
. d r2 1
diS
- ^ s + x l i r r s - L dt-
,
(4-110)
T) = 0,
(4.111)
=
(4 . 1 1 2 )
0 ,7 7
C=
um7)
( •
)
(4.108)
C = 0,17 = 1,
T 2 = T6u),
^-ik(l-S )cf
n
2
2 ) 1^2
7
Ti(77) = T 2(C) = 0,
<* ■
< 1-
0 ,7 7
1,
=
1,
(4.113)
C=l,
(4-114)
^
1 1 P v
( - o - i - 1-
<4-115>
T his is th e full nondimensional system to be analyzed.
4.5.1
N um erical Im plem entation
T h e im portance of making th e transform ation to fix the boundaries is th a t a finite
difference scheme can be employed to solve th e system w ithout the need to ad ap t the
mesh in space to accom m odate th e moving boundary. A m ain objective in this work
is to track th e m elting front S w ith tim e. Since th e problem is one dimensional, the
trajecto ry of the front is expected to be sm ooth w ith time. Im plicit finite difference
schemes are generally efficient methods to track th is front. Note th a t there are two
boundary value problems to be solved; one for the liquid region and one for th e solid
region. We employ a C rank Nicolsen im plicit scheme for each region w ith explicit
reference given to th e moving front S. The intervals
0
<
77
<
1
and 0 < C < 1 are
discretized into N subintervals or IV + 1 nodal points where i denotes the index to
each of these points. Here A z = 1 / N denotes the spatial step size and A t denotes
th e tem peral step size, and
j
denotes th e index in tim e
t.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
70
4.5.2
T he Case o f C onstant D ielectric P rop erties (G eneral Stefan
N um ber)
W ith th e assum ption of constant dielectric attributes, e\, e2, e",and e'2 are unity.
Observe th a t equations (4.104) and (4.105) are decoupled from th e heat equations
in the solid and liquid regions. Using these, together w ith th e boundary conditions
(4.111-4.114), we can deduce an exact solution for th e electric fields in b oth the liquid
and solid portions of the m aterial. The solutions are given as:
(4.116)
0i = A eiksriT>+ B e ~iksriv,
02
= C e ^ 1- 5^
+ D e -W -V ™ ,
(4.117)
where
T] = y/a i + icri
1^2 — n/£c*1 + iSox
A and C are determ ined from the following system:
where
Fix = eikSTl + I i ± i e _ifcsri
p
2i
—piksri _ ri+i -ifcsri
e
ri-ie
b = 2 e~iksr' f ^ - [
W ith analytical solutions of th e electric fields given, the semi- im plicit Crank Nicolson
discretizations need to be carried out only on th e energy equations in each of the
respective (liquid/solid) regions. The discretization of th e system (4.106)-(4.107) is:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where the discrete representations for the spatial derivatives and th e time derivative
are:
df _
dz ~
2Az
d2f
dz2—
df _
dt
f i - l —2 f i + f i + l
(A z)2
At
The coefficients are given as:
*
=
=
C. = - M S f M ~
^
-(1
(M
D, = 4
= (1 - S ) f (C - 1 ) £ -
C2 =
4
2
- S ) f (C - l ) g -
(4.118) and (4.119) hold for
^
2^
2
* 1+
;$ > )
- i^ ,)
B2 = 4
((1
- S )2 +
A
((1
- 5 )2 - , $ , )
= 4
< i < N . For i — 1, boundary conditions (4.108) and
(4.109) yield the following discrete equations:
T f t ' i B i - 2AzSf3Al) +
+ C X) =
+ 2 A z S 0 A x)
~ T \fi{A \ + C i) + P|<£i|2 ,(4.120)
T2,i =
0
(4.121)
The centered differenceformula for the first derivative introduces a fictitious point
T/o
a t the boundary rj =
0.
This point can be represented in term s of T ( i and T (
2
through th e use of the mixed boundary condition (4.108). Observe th a t the centered
difference representation for first derivatives has 0 ( ( A z ) 2) accuracy whereas the one­
sided representation has O ( A z ) accuracy. For i = N + l , boundary conditions (4.109)
and (4.110) yield:
T i ,N + 1 = 0 ,
(4.122)
T { N^ = T bw,
(4.123)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
72
Equations (4.118-4.123) represent the closed discrete system to be solved for both
regions (liquid and solid). In each region, th e discrete system is set up as a tridiagonal
m atrix system to be solved.
4 .5 .3
T h e C a se o f T e m p e r a tu r e D e p e n d e n t D ie le c tric P r o p e r ti e s (G e n e r a l
S te fa n N u m b e r)
In general, the dielectric attributes of a m aterial do depend on tem perature. In this
case the system (4.104-4.115) becomes more complicated to solve as equations (4.104)
and (4.105) governing the electric fields in the liquid and solid regions respectively
are now coupled to the energy equations (4.106) and (4.107). We m ust therefore
resort to numerical methods to solve for these fields along w ith the tem perature
distributions. We use Richmeyer’s M ethod (see Smith G. (1993)) in dealing with
the nonlinear terms involving the tem perature distributions of the liquid and solid
regions. We now discuss the philosophy of the method.
In discretizing the system (4.104-4.115), we use the same finite difference repre­
sentations for the space and time derivatives as before. Notice th a t (4.106) and
(4.107) contain nonlinear term s in 7 \ and
7-2
respectively. In general, we can expand
a nonlinear term (say F ( T) ) as F ( T ) = F ( T j ) -I- ^ ^ ( T - T J) + 0 ( A T 2). This is
ju st the Taylor series representation of F about TJ. Since A T = T J + 1 —T J is small
for small tim e steps, the Taylor series expansion is a good approxim ation. We can
now rew rite F( T ) as F ( T j ) + dFj ^ (T J+ 1 —T j ). Using this finite difference repre­
sentation for the nonlinear term s of (4.106-4.107) along w ith th e finite difference
representations given previously yields the following finite difference scheme:
4>i,i(Ei -F 2zAz/;<S) -I- 2<^ii2 —AikS,
(4.124)
(4.125)
(4.126)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
73
(4.127)
<t>l,N+l ~ 4>l,N = 4*2,2 ~ 4*2,1,
02,1—1 + 4*2^2
4*2 ,n + i (E 2
0 2 ,i- t- 1
+
2 < i < N,
=0,
(4.128)
— 2 i A z k ( l — S )) + 2 ^ 2 ,n = 0,
T f f i B r + Fltl - 2A z S 0 A 1) + T { ^ l ( A 1 + Cx) = ^ ( A
(4.129)
+ Fu + 2 A zS 0 A l)
-TU Ai+ C x)
+ P e " ( T / i ) | 0 i j l |2,
(4.130)
T££ xA, + I ff'(ft + Fw) +3SJ.C. = -Th-iAi + T t/D t+ P u )
- T i ,M C , + P e’U T ’J f a t f ,
2 <i<N,
(4.131)
T i,*+i = 0,
(4.132)
T2A = 0,
(4.133)
T it-,A 2 + T £ \ B 2 + F2i1) + T i ^ C 2 =
- T i l_lA2 + T il{D2 + F2J
- T i i+lC 2 + P 6^{T i{)\4 > 2 ,i\\
2 <i<N,
T2)n + i = Tin,,,
(4.134)
(4.135)
where
E 1 = - 2 + (A z f S W i a ^ T l j + u n e ? ® )
F M = -2AiP^e''(2Y,i)|01,i|2
F2 = - 2 + (A*)2(l - S n H S a ^ i T i j + i S a ^ iT Z J )
Fu = -2AiP^e^(7?)i)|02,i|2
As in the case of constant complex perm ittivity, the problem is posed as tridiagonal
m atrix system
to be solved in the liquid and solid regions of the m aterial. The
tridiagonal entries of the m atrix are different due to
th e nonlinear term s involving
the tem perature field for each region.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
74
035
tim« (I)
F ig u re 4.8 M elting front position S vs. tim e t (constant complex dielectric perm it­
tivity): /3 = 1; k = 1.
Ui
I
s
F ig u re 4.9 M elting front position S vs. tim e t (tem perature dependent complex
dielectric perm ittivity): P=400; /? = 1; k = 1 .
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
75
Figure 4.8 gives the trajectory of the melting front in tim e for P = 400 and
P = 450. Here the complex dielectric perm ittivity is assumed constant. As suggested
by the corresponding steady-state curves, th e melting front for P — 400 propagates
out to S eq = .15. But, for P = 450 the front moves much farther before coming
to rest (S eq = .49). Figure 4.9 gives the trajectory of the melting front in time for
P = 400. Here the tem perature dependence of the complex dielectric perm ittivity is
taken into account. The model incorporating this tem perature dependence predicts
th a t melting front front will propagate much farther to a value S eq = .51 as compared
with the constant perm ittivity model.
Figure 4.10 gives plots of the trajectory of the m elting front position S ( t )
v s.
t
for the case of constant complex dielectric perm ittivity. The solid curve denotes the
trajectory calculated using the leading order asym ptotic result for th e large Stefan
number limit (L —» oo). The dashed curve shows the trajectory for L = 33 and the
dotted curve gives the trajectory for L = 100. These curves are calculated using the
finite difference routine discussed. Note th a t as the Stefan num ber L is increased, the
trajectories given by the finite difference code converge to the asym ptotic result. In
calculating the melting front trajectory, it is far more efficient to employ asym ptotic
approxim ations when L »
1 than it is to use the finite difference algorithm in th a t
the com putational run tim e can be reduced dramatically. Note th a t for L »
1 the
asym ptotic approxim ation serves as a good check on the finite difference algorithm .
Figure 4.11 compares the m elting front trajectory in tim e for a layer of m aterial
whose bottom interface is a dielectric interface, w ith a m aterial whose bottom
interface is prefectly conducting. For th e power setting studied here ( P = 600), the
m elting front propagates faster and farther for the case of a perfectly conducting
interface at (z = 1 ) than it does for th e case of a dielectric interface a t (z = 1 ).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
76
large Stephan num ber limit
U 33
U 100
tima
F ig u re 4 .1 0 M elting front position S vs. tim e
? 0 .4
t
:
P — 600; p = 1; k =
1;
Tbw = —.5.
dielectric interface (z«1)
perfectly conducting interface (z»1)
F ig u re 4.11 M elting front position S vs. tim e t: P = 600; P = 1; k = 1; Tbw = —-5-
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
77
4.6
D iscussion
T he model posed gives some physical insight into th e melting of solids th a t was not
known a priori. For example, th e numerical experiments conducted suggest th a t the
final resting position of the m elting front depends, in a highly nonlinear manner, on
th e power. If it is desirable to m elt only a certain portion of th e solid m aterial, the
model suggests how the power should be controlled in order to achieve the desired
steady-state m elting front position. For example, the steady-state curve (4.7) implies
th a t given an initial melting front position 5j„it = . 1 , in order to move th e front out
to (say S eq = .4), the power term P should increase beyond PCTu = 408, and then
decreased to a value of about P = 300. This is an interesting hysteresis phenomenon.
T h a t is, the trajectory of S eq obtained by increasing P from a value less than Pcru to
a value greater than Pa-it, is not the same as th a t obtained when decreasing P from
a level greater th an Pcru to a level less th an Pa-it- M athematically, when P starts out
below Pa-it, S eq tracks one stable branch of the steady-state curve. Once P exceeds
P = Pcrit it jum ps to another stable branch of this curve, and reducing P forces S eq
to continue along this new stable branch. Note th a t in a transient problem, given
th e initial conditions, the times taken to reach a given state differ depending on the
branch traversed. This provides useful quantitative inform ation for applications.
As was seen in Figure 4.11, th e am ount of power necessary to achieve a desired
m elting front position can be reduced, in many cases, by introducing a reflective
surface a t the bottom of th e m aterial layer.
In commercial processes where it
is desirable to use the least am ount of power possible in processing a m aterial,
th e results discussed may have potential applications. In estim ating power levels
necessary to achieve a desired steady front position, it is im portant, for certain
m aterials, to account for th e tem perature dependence of th e complex dielectric
perm ittivity.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
A P P E N D IX A
THE M E T H O D OF M ULTIPLE SCALES
We present the full equations in dimensional form :
A E ' - V (V • E ') = no
( A l)
P C p [^- + (u; • V)T'] = ktA T ' + | | E
7^7
(yt
+ (u ' • V )u ' =
(A 2 )
'|2
V-— - p /p 0g + i'A u '
Pq
(A 3)
u ' = (u\, u'2, u'3) represents the velocity vector field, T ' is the tem perature P ' is the
pressure and E ' is the electric field,
po, cr and a
are magnetic permeability, the
perm ittivity and the conductivity, respectively, of the dielectric. kt is the therm al
conductivity, p is the density of the dielectric and po is the reference density at the
am bient tem perature, v is the kinematic viscosity of media.
In the above system we recognize th a t there are two tim e scales present in the
physical problem. The oscillatory tim e of the electric field and the diffusive time
for the tem perature. In order th a t we may study the therm o-dynamical aspects of
the physical system we wish to study the long tim e (diffusive time) behavior of the
system.
To do this we scale th e dimensional oscillatory tim e (say ii) w.r.t.
electrom agnetic wave frequency and note th a t the slow time variable (say r ) =
where
8
— k/(ojcL2) is a small param eter. Note ^
the
8 t\
We introduce the
following scalings to nondimensionalize A.1-A.3:
x = x '/d ,
0=
( - 1
t = ait', u = du'/re, P = d 2 / (K2 pa)P '
+ &>
E = B '/E o
R = (agd 3 0o)/{Kv),
P r = is/K,
/3 — h d /k t, x = a 0 \E 0 \2 d2 / ( 2 k td0)-
(A.4)
78
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
79
We arrive a t th e following nondimensional system with the small param eter (6 ):
A E - V (V • E) = k 2
dT
^ + 5{u • V)T] =
du
— +
dt
6
(-4-5)
S ^ +S l (/(r)E)
6 [AT
(-4-6)
+ x / C n |E |2]
(u • V )u = 5 ( - V p — (d 3 / n 2g + R p r T )k -+- pr A u )
(-4.7)
We now perform asym ptotics on the above system to a tta in a leading order
system in 5 th a t is uniformly valid for all time. Let u = u 0 + <5ui, E = E 0 + (5Ei and
T = T 0 + 6 T 1.
The leading order system:
A E „ - V (V • E 0) = k 2
W (‘r K ) +
dT 0
=
dt
duo
=
dt
(-4.8)
^ I (/(T“)E',)
0
(-4.9)
0
(>1. 1 0 )
V -u = 0
(-4.11)
Hence Ta = p(x, r) , u 0 = h (x , r) and E 0 = L ( r ) E 0 (x)eiu'tl are solutions.
1« L = ± § ; U ' ( T°)T iE .) + g r O P V lT .E J
M = e & (e .( X .) E„) + ^ - i ( / ( T „ ) E 0)
The 0 ( 5 ) system:
AE! -
V (V • Ex)
^
at
= ^ [ ^ ( f r E x ) + ^ A ( / ( T 0 )E i) + L + M]
= - ( u c • V)T 0 - ^
= — (Uo • V )u 0
ar
- ^
+ A T 0 + x /( T 0 ) |E 0 |2
(A. 12)
(-4-13)
+ ( - Vp - (d3 / n 2g + RprT)~k + p r A u 0)(A 14)
V ■Ui = 0
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(A 15)
80
Integrating A. 13 w ith respect to the fast tim e t and letting t go to infinity we
get:
- ( u 0 • V )T 0 - ^
OT
+ AT 0 + x f ( T 0) lim i
t-* o o t JO
= 0
(A.16)
note th a t u\ = T(x, r ) t where T denotes the right hand side of A.14. In order for
our leading order expansion to be uniformly valid over all tim e we require th a t T
equal zero.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
A P P E N D IX B
D E R IV A T IO N O F E Q U A T IO N S G O V E R N IN G “ C O M P L E X
E L E C T R IC F IE L D ”
From Maxwell’s equations we deduce the following general equation governing the
behavior of the electric field within in dielectric medium.
A E - V ( V - E ) = - ^ ( er( r ) E ) + | (
c
2 2
)E )
(B .l)
From our two-timing analysis we determined th a t a and er only evolve on the slow
tim e scale r . Hence, we can place these variables out of the tim e derivative operators
to get:
(B.2)
A E -V (V -E ) = -
We expect the electric field to have a time harmonic behavior so th e solution is
assumed of the form:
E (x , t) = A (x)sin(w l) + B(x)cos(cjf)
(B-3)
S ubstituting this expression for E into B.2 and equating sine and cosine term s yields
two equations:
A A - V (V - A) = - ( - e r ( 7 > 2 A ) - u ; ^ ( B )
c
to
(B.4)
A B - V (V • B ) = - ( —er ( 7 V B ) + w ^ ( A )
c
(B.5)
Let e Q = A + i B. Adding i times B.5 to B.4 gives:
A e 0 - V (V • Co) = - k 2 [er (T )e 0 + ze"(T)e0]
81
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(B.6 )
82
where e"(T) =
note the physical electric field E is represented in term s of e Q as
E = real(e 0 e- “‘,t)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
A P P E N D IX C
W K B ANALYSIS
We consider the following ordinary differential equation :
rf- 'Ij
—
dz
~XQ(z)u = 0
(A = —?)
o
u= 0
where Q(z) =
0
<
(C .l)
2 < 1
z = 0,1
(C.2)
+ a2
is a mono tonic decreasing function of z which changes sign in
Note th a t
the interval (0 , 1 ) and th a t a 2 is bounded by
a t .z =
1,
otherwise th e solution
to the above system would be trivial. Since a 2 m ust be positive we have a turning
point for the above equation in the interval given. As A —> oo we can apply WKB
theory to get solutions to
2 .6
away from the turning point.
On the interval 0 < = z < zc where zc is the turning point, u has the form:
u =
Voi^dz) + _ J L _ C( - A r \/5(I)*‘)
Q( z ) t
Q(z)i
(C.3)
To satisfy condition th a t u/(0) = 0 requires th at u / take the form:
Ul =
A
i C(a
l:c Voiijdz) _
e ( - A / ; c y/otr)dz+2xf;c y/odm
(c .4 )
Q(z)*
b ut this solution blows up away from the boundary as A —> oo. Therefore, the
solution u i is equal to zero
Let uj[ denote th e solution ju s t about the turning point. To study behavior
of the solution about turning point let z' = (Q’(zc)) 3 ( 2: — zc)A§. The governing
differential equation for u in term s of z ' becomes:
0 - A - O
83
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(C.5)
84
u // = BAi ( z ' )
(C. 6 )
Note th a t Ai(z' ) —» 0 as 2 ' —►oo which is consistent with our conjecture on uj (see
Bender and Orszag). As 2 ' —» —oo u // looks like:
C s i u ( U - z ’) l + f )
< C ' 7 )
The solution u u i of region III [zc < z < 1) has the form:
uni =
^ r s in { \ [ J - Q ( z ) d z ) + ^
. cos(A [ s j - Q { z ) d z )
- Q ( Z ) ‘i
-< 2(z ) 4
^
(C. 8 )
The boundary condition u m { 1) = 0 requires th a t u u i have the following form:
unr =
^
L sin{X f
—Q ( z )4
We m ust now match ur n as
moves toward
2
2
J-Q{z)dz)
(C.9)
Jz
moves toward the turning p t (z = zc) w ith
uji as 2
= 1. The solution u u i can be w ritten as:
«/// =
c_
~ Ls in { \{ [ \J—Q( z ) dz - f \ f - Q ( z ) d : ) )
Jzc
JZr
Jzc
3 ( 2 ) 54
—-<CJ(Z\
(C.10)
A m atching of the Airy solution u u with u u i requires th a t the constant expression,
^ fzc \ J ~Q{ z ) dz =
2 n 7r -1-
\ for A > > 1. (n = ------- 2, - 1 , 0 , 1 , 2 • • •). Hence, 2 C m ust
approach 1 for large disturbance wave numbers. For zc — z «
1, A /Zc yJ—Q( z ) dz
can be approxim ated and w ritten in term s of z' as 2 /3 (—z ' ) i . So near the turning
point U u exhibits th e correct behavior to m atch w ith th e Airy function. Recalling
th a t 2 C —> 1, the integral / z* yJ—Q( z ) dz can be approxim ated by Taylor expanding
2 a?
Io~2) ^
Q{z) about 2 = 1. The approxim ation in term s of 9 and a i s
----- Recalling
th a t this expression equals a constant, we arrive a t the following equation:
j 2 /3
4
CT3
- 071/ 0 ~ ^
O4/ 3
Note as a —> 0 0 a 2 —>
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
( C . ll )
REFERENCES
1.
C handrasekhar S. 1981 Hydrodynamic and Hydromagnetic Stability. Dover
Publications Inc., New York, NY.
2. Drazin P. G., Reid W. H. 1981 Hydrodynamic Stability. Cambridge University
Press, New York, NY.
3. Fang C. S., Lai P. M. C. 1995 Microwave heating and separation of water-in-oil
emulsions. Journal o f Microwave Power and Electromagnetic Energy. 30,
46-57.
4. Hippel A. V. 1954 Dielectric Materials and Applications. John W iley and Sons
Inc., New York, NY.
5. Kriegsmann G. A. 1992 Therm al runaway in microwave heated ceramics: A
one-dimensional model. Journal o f Applied Physics. 71, 1960-1966.
6
. Kriegsmann G. A. 1993 Microwave heating of dispersive media. S IA M Journal
o f Applied Mathematics. 53, 655-669.
7. Kulacki F. A., Goldstein R. J. 1972 Thermal convection in a horizontal fluid layer
with uniform volumetric Energy Sources. Journal o f Fluid Mechanics. 55,
271-257.
8
. Pellow A., Southwell R. V. 1940 On convective motion in a fluid heated from
below. Proc. R. Soc. A 176, 312-343.
9. Ramo S., W hinnery J. R., Van Duzer T. 1984 Fields and Waves in Com m uni­
cation Electronics. John Wiley and Sons Inc., New York, NY.
10. Eds. Donnelly, R. J., Herman, R., Prigogine, I. 1966 Non-equilibrium Thermody­
namics, Variational Techniques and Stability. University of Chicago Press,
Chicago, IL.
11.
Roberts, P. H. 1967 Convection in horizontal layers with internal heat generation
theory. Journal o f Fluid Mechanics. 30, 33-49.
12. Sparrow, E. M., Goldstein, R. J. and Jonnson, V. K. 1964 Therm al insta­
bility in a horizontal layer: Effect of boundary conditions and non-linear
tem perature profile. Journal of Fluid Mechanics. 18, 513-528.
13. Stuerga D., Lallem ant M. 1993a An original way to select and control
hydrodynam ic instabilities: Microwave heating. P art I: Hydrodynamic
background and th e experimental set-up. Journal o f Microwave Power
and Electromagnetic Energy. 28, 206-218.
85
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
86
14. Stuerga D., Lallemant M. 1993b An original way to select and control
hydrodynamic instabilities: Microwave heating. P art II: Hydrodynamic
behavior of water and ethanol under microwave heating and reduced
pressure. Journal o f Microwave Power and Electromagnetic Energy. 28,
219-233.
15. Stuerga D., Lallemant M. 1994 An original way to select and control hydrodynam ic instabilities: Microwave heating (P art IV). Journal o f Microwave
Power and Electromagnetic Energy. 29, 3-30.
16. Stuerga D., Lallemant M., Steichen-Sanfeld A. 1994 An original way to select
and control hydrodynamic instabilities: Microwave heating (P art III).
Journal o f Microwave Power and Electromagnetic Energy. 29, 3-19.
17. Stuerga D., Lallemant M. 1994 An original way to select and control hydrodynam ic instabilities: Microwave heating (P art IV). Journal o f Microwave
Power and Electromagnetic Energy. 29, 20-30.
18. Tveitereid, M. and Palm , E. 1976 Convection due to internal heat sources.
Journal of Fluid Mechanics. 76, 481-499.
19. W atson P. M. 1968 Classical cellular convection w ith a spatial heat source.
Journal of Fluid Mechanics. 32, 399-411.
20. Yucel A., Bayazitoglu, Y. 1979 Onset of convection in fluid layers w ith nonuniform volumetric energy sources. Transactions o f the ASM E, Journal
o f Heat Transfer. 101, 666-671.
21. Hill Jam es M., Dewynne Jeffrey N. 1987 Heat Conduction. Blackwell Scientific
Publications., Boston, MA.
22. Crank J. 1984 Free and Moving Boundary Problems. Clarendon Press., Oxford,
England.
23. Stuerga D., Zahreddine I., More C., Lallemant M. 1993 Bistable behavior in
microwave heating: T he first experimental evidence. C. R. Acad. Sci.
Paris, t. 316, S eries II , 901-906.
24. Nachman M., Turgeon G. 1984 Heating p attern in a multi-layered m aterial
exposed to microwaves. IE E E Transactions on Microwave Theory and
Techniques, M T T -3 2 N O . 5, 547-552.
25. Streit U. 1991 One-dimensional implicit moving boundary value problems: A
new front tracking approach. Numer. F und. Anal, and Optimiz., 12 (2
a n d 4 ), 413-428.
26. Sethian J. A., Strain J. 1992 Crystal growth and dendritic solidification. Journal
o f Computational Physics., 98, 231-253.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
87
27. Sm ith G. D. 1993 Numerical Solution o f Partial Differential Equations: Finite
Difference Methods. Oxford University Press Inc., New York, NY.
28. Burden R. L., Faires J. D. 1989 Numerical Analysis. PW S-KENT Publishing
Company, Boston, MA.
29. R ubinstein L. I. 1971 The Stefan Problem. American M athem atical Society,
Providence, RI.
30. G oodm an T. R. 1958 The heat-balance integral and its application to problems
involving a change of phase. Transactions o f the ASM E., F e b ru a ry , 335342.
31. Boley B. A. 1961 A method of heat conduction analysis of melting and solidifi­
cation problems. J. Math. Phys.., 40, 300-313.
32. C rank I. 1957 Two methods for the numerical solution of moving-boundary
problems in diffusion and heat flow. Quart. J. Mech. Appl. Math.., 10,
220-231.
33. Douglas J., Jr. 1957 A uniqueness theorem for the solution of a stefan problem.
Proc. Am er. Math. Soc.., 8 , 402-408.
34. Friedm an A. 1959 Free boundary problems for parabolic equations. I, Melting
of solids. J. Math. Mech., 8 , 499-517.
35. Feynm an R.P., Leighton R.B., Sands M. 1964 The Feynman Lectures on Physics.
Addison-Wesley Publishing Company, Reading, MA.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
IMAGE EVALUATION
TEST TARGET ( Q A - 3 )
im
1.0
|06
■AO
m
12.0
l.l
1.8
1.25
1.4
1.6
150mm
IIVMGE . I n c
1653 East Main Street
Rochester, NY 14609 USA
Phone: 716/482-0300
Fax: 716/288-5989
© 1993, Applied Image, Inc., All Rights Reserved
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Документ
Категория
Без категории
Просмотров
0
Размер файла
3 584 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа