# 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.

1/--страниц