r a a 1i t B National Library or Canada Canadian T h e s e s S ervice Biblioih&que nationale du Canada S ervice d e s th e s e s can ad ien n es Ottawa, Canada K1A 0N4 NOTICE AVIS The quality of this microform is heavily dependent upon the quality of the original th esis submitted for microfilming. Every effort h as been m ad e to en su re the highest quality of reproduction possible. La quality de cette microforme depend grandement de la quality de la t l i s e so u m ise au microfilmage. Nous avons tout fait pour assurer une qualit6 sup6rieure de reproduc tion. If pages are missing, contact th e university which granted the degree. S'il m anque des p a g es, veuillez communiquer avec I'universit6 qui a conf6r6 le grade. S om e p a g e s may h ave indistinct print especially if the original p a g e s were typed with a poor typewriter ribbon or if the university sent u s an inferior photocopy. La quality d'impression de certaines p a g es peut laisser a d^sirer, surtout si !es p a g e s originales ont 6t6 dactylogra p h i e s k I'aide d'un ruban u s6 ou si I’universit6 n ou s a fait parvenirune photocopie d e quality in irieu re. Reproduction in full or in part of this microform is governed by the Canadian Copyright Act, R.S.C. 1970, c. C-30, and subsequent am endm ents. La reproduction, m em e partielle, de cette microforme est soum ise k la Loi ca n a d ien n e sur le droit d'auteur, SRC 1970, c. C-30, et s e s am endem ents subsSquents. NL-339 It. 88/04) c R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. Canada W IR E GRID M ODELING FOR M ICRO W A V E HEATING AND THERM AL RUNAW AY by Tien D. Pham A thesis submitted in conformity with the requirements for the Degree of Master of Applied Science in the University of Toronto Department of Electrical Engineering University of Toronto Toronto, Ontario © Tien D. Pham 1990 R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. eg H S . B National Litxary of Canada Bibliotheque nationals du Canada Canadian Theses Service Service des theses canadienoes O ttaw a. Canada XIA ON 4 The author has granted an irrevocable non exclusive licence allowing the National Library of Canada to reproduce, loan, distrfoute or sell cop ies of his/her thesis by any means and in any form or format, making this thesis available to interested persons.' L’auteur a accorde une licence irrevocable et non exclusive permettant & la Bibfiothdque nationale du Canada de reproduire, prSter. distribuerou vendre des copies d e sa th ese de quekjue manidre et sous quelque forme que c e soit pour mettre d es exemptaires d e cette these a la disposition des personnes interessees. The author retains ownership of the copyright in his/her thesis. Neither the thesis nor substantial extracts from it may be printed or otherwise reproduced without his/her per mission. L’auteur conserve la propriety du droit d’auteur "qui protege sa these. Ni la these ni d es extraits substan tiets de celle-ci ne doivent etre imprimes ou autrement reproduits sans son autorisation. ISBf! 0 - 3 1 5 -6 5 4 4 8 -1 R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. THE UNrVERSTTY OF TORONTO LIBRARY MANUSCRIPT THESIS AUTHORITY TO DISTRIBUTE NOTE: T h e AUTHOR will sig n in o n e o f the tw o places indicated. It is the intention of th e University that there b e NO RESTRICTION on the distribution of the publication of t h e s e s s a v e in exceptional ca ses. (a) Im m ediate publication in microform by th e National Library is authorized. Author’s signature ____________________ ^ C Date--------- ------------— lA — (JLj-ii;-- -O R - (b) Publication by the National Library is to be postponed until_______________________ (normal maximum d ela y is tw o years). Author's sig n a tu r e __________________________________________ 19 D ate This restriction is authorized for r e a s o n s which seem to m e, as Chairman of th e G raduate D epartm ent of ______________________________________________ , to be sufficient. Signature of G raduate D epartm ent C h a irm a n __________________________________________________________ D a t e ________________________________________ B O R R O W E R S undertake to give proper credit for any u se m ade of the th e s is , and to obtain the co n se n t of the author if it is proposed to m a k e e x te n siv e quotations, or to reproduce th e th e sis in w hole or in part. A ddress Signature of borrower I I i I - I I I I I R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. Date | ABSTRACT The wire grid modeling method based on the thin wire program by J.H. Richmond has been modified to model both conducting and dielectric bodies. The method is vali dated bv comparing the numerical and analytical solutions for the field and power distri butions in a dielectric-filled terminated waveguide model. The method is then used to invest gate the microwave heating of some dielectric materials under various cir cumstances. Computation of cower dissipation in thin dielectric bodies with high loss factors simulates some commonly observed phenomena in microwave heating like edge heating and thermal runaway, which can be satisfactorily explained by simple surfacewave theory. One-dimensional microwave heating in a half-space has also been studied for several dielectric materials to give some insights into the thermal runaway effect. A numerical method is utilized to solve the coupled wave-thermal equation, which appears in the problem as a result of the temperature dependence of the dielectric constant and loss factor o f many dielectric materials. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. ACKNOWLEDGMENTS I wish to express my sincere gratitude to my supervisor Professor Keith G. Balmain for his expert guidance and enthusiastic support during the course of this work. I must also thank Dr. Colin C. Bantin and Dr. Mark A. Tilston for many stimulating discussions and suggestions on antennas and electromagnetic theory. The assistance of Dr. Claude Lorenson and Christine Gallemeault of Alcan International Ltd. is also greatly appreci ated. I would like to thank everyone in the lab who has helped me in many ways. They have provided the atmosphere and knowledge that have made studies at this university an enjoyable experience. Finally, I would like to thank Professor Balmain and Alcan International Ltd. for providing financial support and technical materials that have made this research project possible. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission. TABLE OF CONTENTS INTRODUCTION................................................................................................................ 1 1. FUNDAMENTAL PR IN C IPLES................................................................................... 3 1.1. Basic theory of microwave heatin g............................................................. 3 1.1.1. Dielectric lo s s ................................................................................... 3 1.1.2. Power dissipation............................................................................. 5 1.1.3. Dielectric properties.................................................. 6 1.2. Thin-wire implementation of the moment method ................................... 9 1.3. Fields in rectangular waveguides and co n tain ers...................................... 10 1.3.1. Fields in rectangular w aveguides.................................................... 10 1.3.2. Fields in rectangular containers ...................................................... 14 1.3.3. Excitation of TE w av es.................................................................... 16 2. WIRE GRID MODELING OF CONDUCTING AND DIELECTRIC BODIES ......................................................................................... 19 2.1. Modeling of continuous media using wire g r id s ........................................ 19 2.2. Dielectric heating in a small, deep rectangular co n tain er......................... 26 2.3. Dielectric heating in other rectangular containers .................................... 34 3. ONE-DIMENSIONAL MICROWAVE HEATING THERMAL R U N A W A Y ....................................................................................... 38 3.1. Wave propagation in inhomogeneous media ............................................. 38 3.2. Numerical method for coupled wave-thermal equation............................ 41 3.2.1. Boundary condition considerations................................................. 42 3.2.2. Numerical approach ......................................................................... 44 3.3. Thermal runaway in one-dimensional microwave h eatin g ....................... 48 3.3.1. Example 1: Cooked beef ju ic e ........................................................ 49 3.3.2. Example 2: Nylon 66 ....................................................................... 51 3.3.3. Example 3: Thawing of frozen water ............................................ 53 R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 4. TRANSVERSE MICROWAVE HEATING - THERMAL RUNAWAY 59 4.1. Surface wave theory ...................................................................................... 59 4.1.1. TE w aves........................................................................................... 61 4.1.2. TM w aves.......................................................................................... 63 4.2. Edge heating in infinitely long dielectric strips.......................................... 64 4.2.1. The TE case ................................................................................. —- 64 4.2.2. The TM case ..................................................................................... 67 4.3. Edge heating in rectangular dielectric s la b s ................................................ 4.4. Thermal runaway in transverse microwave heating................................... 68 73 CONCLUSIONS............................................................................................................. 78 .APPENDICES A IMPLEMENTATION OF ELECTRIC AND MAGNETIC SYMMETRIES IN THE THIN-WIRE PROGRAM ....................................... 79 B COMBINATION OF SYM M ETRIES..................................................................... 83 REFERENCES ................................................................................................................ 86 R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission. 1 INTRODUCTION A widely-used application of the electromagnetic moment method is the wire grid modeling of two-dimensional or three-dimensional bodies in electromagnetic scattering problems. By this method, a solid body is simulated by a grid of wires, from which many electromagnetic characteristics of the body can be determined [1], [2]. The method has been widely used to compute, for example, radar cross sections of aircraft and fields of parabolic reflectors. In general, the main use of the wire grid method has been the model ing o f metallic-surface bodies, followed by computations of far field radiated from the bodies. The application of this method can be extended significantly if one ca^ model bodies of other materials, rather than only perfect conductors. The motivation for this idea has its origin in the electromagnetic problems of microwave heating, where one wants to determine the power dissipation in food when exposed to microwave energy. Basically, microwave heating problems are similar to scattering problems involving con ducting bodies. There are, however, several differences between these two classes of problems. First, in microwave heating problems, one is interested in finding the near field or the induced current that gives rise to power dissipation in a lossy dielectric body. Usu ally, near field computation is more demanding than far field computation with respect to numerical accuracy [3]. Second, scattering problems involving conducting bodies are essentially simpler than microwave heating problems, because the zero tangential electric field boundary condition can be applied on the surface of a conductor in most practical cases. The fact that no comparably simple boundary condition exists on the surface of a dielectric body makes the field computation inside it a relatively difficult task. Neverthe less, it will be shown that a dielectric body can be successfully modeled by a wire grid properly loaded with lumped impedances; furthermore, the wire grid modeling method offers the unique flexibility of using wire-, surface- and volume elements together to sup port the modeling of most practical conducting and dielectric bodies. Therefore the R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 2 method is suitable for the analysis of microwave heating problems, including uneven heating and thermal runaway in food materials, which are commonly observed in domes tic and commercial microwave ovens. Since microwave heating is a relatively new and unexplored research area, some basic investigations on the runaway effect will be first conducted for the one dimensional case with the hope that this will provide some insight into the phenomenon. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission . 3 1. FUNDAMENTAL PRINCIPLES This chapter contains a review of the theory of microwave heating along with a short summary of the thin-wire implementation of the moment method. Also, the theory of a waveguide filled with lossy dielectric will be introduced. These fundamental princi ples serve as a basis for the discussions taking place in the chapters that follow. 1.1. Basic theory of microwave heating It has long been known that an insulating material can be heated by applying energy to it in the form of electromagnetic waves. The origin of this heating lies in the ability of the electric field to polarize the charges in the material and the inability o f this polarization to follow the rapid reversals of the electric field. Coupled w; h the polariza tion effects, a dielectnc material can be heated through direct conduction effects due to, for example, the redistribution of charged particles under the influence o f the externally applied electric field. The theory of such polarization and conduction effects has been extensively studied and summarized in the past [4]. In this chapter, however, only a short and simple formulation will be presented to give the reader a basic understanding of the theory of microwave heating. 1.1.1 Dielectric loss At microwave frequencies (300 MHz < / < 30 GHz), dielectric loss in a food material is due mainly to dipolar polarization in which the permanent dipoles contained in the material reorient under the influence of a changing electric field. Another source of polarization arises from charge build-up in interfaces between components in heterogene ous systems, termed Maxwell-W agner polarization. Other kinds o f polarization are elec tronic and atomic polarizations, which dominate at v e r v high frequency (optical) and are not involved in microwave heating. For all kinds of polarization, dipole charges (bound R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission. 4 (bound charges) are induced and the field vector that relates to the bound charges is termed the polarization P. The electric fiux (displacement) density vector D is given by D = EoE + P (1-1) D = £r EoE (1-2) P = (er - l ) e o E (1-3) By definition, which together with eq. (1) yields In general, cr is a complex constant due to the phase difference between P and E. Furth ermore, the polarization vector P lags the electric field vector E, so the relative complex permittivity must be expressed by £. = er - j e p (1-4) where £r and £p are both positive and real. It should be noted that no conduction effect has been taken into account so far. For a dielectric having the conductivity a , the total current density is given by M axwell’s equation VxH = J = o E + - ^ dr where c E isthe conduction current density and 3D (1-5) isthe displacement current density. For sinusoidal electric field variations with an em time variation, eq. (5) can be rewrit ten as J = a E + y'o)£o£rE = ;cd£o (e'r - j Z p ’ - } - £ - ) E GjEO (1-6) M . . . In eq. (6), a/co£o represents conductivity loss and ep polarization loss.Since, with most dielectric measuring techniques, it is difficult to separate the losses due toconduction and R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 5 those due to polarization, it is most convenient to group them together as the loss factor which, in general, replaces ep in eq. (4). The ratio of the loss factor to the dielectric constant is the loss tangent ( 1-8) 1.1.2 Power dissipation Microwave heating involves the conversion of electromagnetic energy into heat. Energy is transported through space or any medium, by means of electromagnetic waves. The power flow through a closed surface S can be calculated by the complex Poynting theorem s V (1-9) The surface integral on the left hand side represents the rate of flow of energy inward through the surface enclosing the volume, the first term in the right hand side represents the stored magnetic and electric energy in the volume and the last term the average dissi pated power. The dissipated power density is given by p = ^coeoei-'E’ -E = y taeo e" | E | 2 (1 10) or in terms of current density J given in eq. (6) ( 1- 11) R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 6 1.1.3 Dielectric properties AT The power dissipated in a dielectric body depends on the loss factor er as well as the electric field E, which is indirectly dependent on the complex permittivity £r . Knowledge of the dielectric properties of the material to be processed by microwave heating is essendal for prediction o f the heating pattern inside the body when it is exposed to microwave illumination. It has long been known that both the dielectric con stant er and the loss factor e r are frequency and temperature dependent A considerable amount of data for inorganic, organic and biological materials has been tabulated on the variation of £r with / and T. For dielectric materials with polar structures, the dependence of £r on frequency / can be explained from the dielectric relaxation of the material. Debye first advanced a theory which expresses £r and Er as functions o f / i n closed forms ' ** due to a given kind of polarization [4], The variation of £r and £r w ith /in the entire fre quency range is then a combined effect of all kinds of polarization with each of them being dominant in a particular frequency band. In the context o f this thesis, we are mainly interested in the dependence of Er and £r on temperature because it has a direct influence on the dynamic heating process and sometimes can lead to the thermal runaway phenomenon. Generally, the chemical bonds in a dielectric material are dependent on temperature and hence change the polarizability and the conductivity o f the material which in turn change the dielectric constant and the loss factor. O f special interest are food materials that usually contain a significant amount o f water. Liquid water is strongly polar in its structure, and therefore able readily to absorb microwave energy and convert it to heat. Therefore water absorbed in food materials has a marked effect upon its dielec tric properties. Bengtsson and Risman [5] showed that there are large differences in £r and £r between ice and liquid water and the dielectric properties o f most foodstuffs fol low a similar pattern to those of water in bulk when plotted as a function o f temperature, as shown Fig. 1-1. The sharp increase of £r around 0°C suggests a tendency to thermal runaway when the food materials are defrosted. As soon as a few droplets of liquid water are produced, the microwave energy is preferentially dissipated in the liquid phase R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission. 7 causing uneven thawing with damaging effect. The temperature dependence of dielectric properties of some typical food materials above freezing point [6] is shown in Fig. 1-2. Because the power dissipated in the dielectric materials also depends on the magnitude of the electric field, thermal runaway can happen in many cases, where the loss factor remains constant but the field distribution is shifted toward the regions of relatively high temperature due to the inhomogeneity in the material. Although the dielectric properties of a material are definitely o f great importance, other factors such as specific heat, ther mal conductivity, mass transfer and especially, the latent heat associated with the phase changes in a thawing process may also affect the heating process significantly. Due to the quick heating process in most cases where there is no phase change, it is practically enough to consider the effects caused by the behaviour of the dielectric properties only. D istilled w a te r 80 iked 7Q [ M ash ed p o ta to 60 I I w 1 50 R aw b eet A 40 ( 30 Cooked 20 C ooked ham \ ^ 10 , y ,.— -20 ______ 25 - 20 T CC) 20 . Cooked v %j C ooked / i ham / . 20 -10 0 20 T CC ) 60 Fig. 1-1 Dielectric properties of some food materials around freezing point [5] R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. 8 • Cooked bee! a Raw ♦ Cooked b e ef 3 0 0 MHz 9 1 5 MHz 2 4 5 0 MHz 15 0 80 juice beef 60 100 u r *— ■ 50 .oA —4—-4*~ T 25 45 T ( °C ) 65 • 5 C o o k e d t u r k e y j ui c e Raw turkey Cooked turkey 25 T 45 ( 5C) 65 * > .' • 100 5 25 45 T (°C ) 65 45 T CC ) 65 Fig. 1-2 Dielectric properties of some food materials above freezing point [6] R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 9 1.2 Thin-wire implementation of the moment method The numerical solutions to electromagnetic scattering problems can be formulated in terms of currents induced in the scattering bodies. A solution for these currents is first sought and then the scattered field is calculated from these induced currents. There are many methods of solving for the induced currents and some of these methods have been summarized by seven 1 authors [7], [8]. Generally, these methods are based on formulat ing an integral equv ;on • : r th'* induced currents. By expanding the current distribution as the weighted sum of the . ter; hers o' some basic set, the integral equation can be used to create a system jF lin • .’ quatior.s from which the current distribution can be solved algebraically. This method o f solution is generally known as the moment method. A special implementation of this method using thin wires was developed and pro grammed by J.H. Richmond [2] to analyze antennas and structures of interconnected wire segments. One of the motivations for the development o f the program was wire grid modeling of conr.mous surfaces, from which the scattering properties of the surfaces can be determined e.vily. Such questions as "How dense should the wire grid be to properly model a given surface ?" or "what is the best wire radius to grid spacing ratio ?" arise automatically in c:-'section with the modeling of wire grids and have been investigated by several authors [9], [10], [11]. As a matter of fact, the program’s effectiveness has led to the publication of many papers concerning the modeling aspects of scatterers. The use fulness of the program is due mainly to the straightforward approach it offers the user when modeling scatterers with bodies of arbitrary shape. In effect, the physical boundary-value problem relating to the electromagnetic scattering problem is reduced to a network problem in which the network’s parameters must be found as the impedance matrix of the scatterer in order to solve for the induced currents. It must be realized that the range of problems solvable on any real computer extends to bodies not exceeding several wavelengths in size. This limitation is due to both the cost of computing time and the limited memory capacity of computers. In an effort to reduce the number o f unknown R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 10 currents in a wire grid model, electric and magnetic symmetries have been utilized in the program. The concept of ehctxic and magnetic symmetries is best illustrated in Jordan and Balmain [12]. Hovvr er, a short treatment of these symmetries together with some practical aspects of t'tei use ivt the program can be found in Appendices A and B. As a step in the evolution o t h e program and its use, volume modeling of three dimensional dielectric bodies will be presented in Chapter 2. The use of this new approach in the rest of the thesis will be concentrated on the computation of power dissipation or near field in various wire grid models. As a final note, the original version of the thin-wtire program by Richmond has been reformulated by M. Tilston [13] to suppress the computation of nonphysical asym metric fields and to allow analysis of multiradius thin wire structures. 1.3 Fields in rectangular waveguides and containers It is seen from eq. (10) that the power density is proportional to the loss factor and the squared magnitude of the electric field. A knowledge o f the field distribution inside the dielectric body is therefore essential in determining the microwave heating pattern. An understanding of the basic principles o f microwave heating can be gained by studying the electromagnetic field distribution inside a rectangular metallic container filled with lossy dielectric. This configuration offers the unique advantage o f having a simple geometry and boundary conditions. In fact, the well established theory of rectangular waveguides can be applied directly in some cases to give the solution for the electric field distribution inside the container. 1.3.1 Fields in rectangular waveguides. The theory of rectangular waveguide can be found in any standard text about elec tromagnetics [12], so only a short summary of the theory will be presented here to facili tate the discussions that follow. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 11 In general, Maxwell’s curl equations and the wave equations are VxH = ;coeE V2E = y 2E VxE = ~ M i H (1-12) V2H = t2H (1 -1 3 ) e = £o (e'r-jE r) (1-14) where y = V -r to2i j r Consider the waveguide filled with lossy dielectric shown in tions in the z direction may be expressed Fig. 1-3.Assume that varia as e -Y' z, where ingeneral j 2= oq + j $ z is the propagation constant in the z direction, yet to be determined. M axwell’s curl equations written in rectangular co-ordinates for the dielectric region within the guide become x Fig 1-3 Geometry of perfectly conducting rectangular waveguide R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. and the wave equations for Ez and Hz are + 72z E z = - co2)i££7 ox2 (1-16) dy2 (1-17) Equation (15) can be rewritten to express Hx, Hy , Ex and Ey in terms of Ez and Hz alone. There must exist either Ez or Hz so that other components do not all vanish. It is con venient to divide the possible field configurations within the guide into two sets, transverse magnetic (TM) waves for which Hz = 0, and transverse electric (TE) waves for which Ez = 0. Assuming perfect conductivity for the walls of the guide, one can set the tangential electric fields at the walls equal to zero: Ex = Ez = 0 at y = 0 and y = b Ey = Ez - 0 at x = 0 and x - a Using the separation of variables technique and applying the specified boundary condi tions to solve for Ez and Hz from eqs. (15) and (16), we get the complete solution for TM and TE waves as follows: R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission. 13 TM waves. —y2E2oa _yz Ex = ------ =— cosAx siaBy e h - y zEz0B _jyZ Eyy = -------= h i --- sinAx cosBy e Ez - E z o sinAx sinBy e Y'z j(3iCEzQB Hx = -----h2 sinAc cosBy e _v , (1-18) -jox.Ez qA . -77 ____ p in //„y = --------= ----- cosAxv sinBy e" l2 h: Hr= 0 n •1 TE waves. jw\iHzQB . _TjZ = -------=----- cosAx sinBy e h2 „ -J' w H zoA . , n _v £ yv = ---------= h2 ----- sinAx cosBy e Ez = 0 Tz//zoA . -yz / / , = -----=— sinAx cosBy e “ h2 yzHznB _YZ / / v = ----- ;— cosAx sinBy e (1-19) , , Hz = HzQ cosAx cosBy e where A =^ ~ a £ = — b m,n = 0,1,2,... and /n+rz > 0 yz = ^ A 2 + B l + f R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. ( 1-20 ) ( 1- 2 1 ) 14 (1-22) h 1 =7^ + C02p£ For waveguides filled with lossless dielectric, yz is either imaginary or real, which corresponds to propagating or cutoff (attenuating) modes, respectively. For waveguides filled with lossy dielectric, yz is always a complex constant, which means every mode propagates and attenuates simultaneously and there is no sharp division between pro pagating and cutoff modes. 1.3.2 Fields in rectangular containers A rectangular container can be considered as a rectangular waveguide terminated at z = l as shown in Fig. 1-4. Due to reflection at the bottom wall, a TM/TE wave travelling in the positive z direction will result in a corresponding wave travelling backward in the negative z direction. The field distribution within the container can be found as sum of these two waves. Recall that in the derivation of TM/TE waves above we have made the assumption that the z-variation of the field within the guide can be expressed as e " , so z x a Fig. 1-4 Geometry of rectangular container: side and bottom walls are perfectly conducting, while the top (z = 0) is open. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 15 that backward travelling waves can be represented by assuming a variation o f the form e lz~. It is easy to see that the transverse (x- and y-) variations must be identical for for ward and backward travelling waves, at least to satisfy the zero tangential electric field condition at z = / for all values of x and y (0 <x <a, 0 <y< b). With this in mind, the zvariation of every E or H component can be expressed as Z(z) = e~y' z + T e y’z (1-23) The reflection coefficient T can be found by applying the boundary condition Ex = 0 at z = /. Then we have Z ( l)^ -e ~ y‘l + T e y‘l = 0 which gives T = -e ~ lyi1 (1-24) Thus, thefield distribution within the container due to a particular mode is the same as given in eq. (18) or (19) with the z-variation e~y,z beingreplaced by Z(z) given in eqs. (23) and (24). The total field is then a combination of all the modes existing within the container. When the container is illuminated by some external electromagnetic source, different modes will be excited depending on the geometry of the container, the material filling the container and the illumination of the container. Determination o f the magni tude for each mode would be a complicated task. It is well known from the theory of waveguides that only a few lowest order modes can propagate in a given waveguide due to the cutoff state o f the higher order modes. In a container, which can be regarded as a terminated and usually short waveguide, the evanescent field due to the cutoff modes can be significant at the opening of the container and thus cannot be ignored. It turns out, however, that for containers with small cross-sectional dimensions, a properly chosen excitation source actually can generate only the lowest (dominant) mode. R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. 16 1.3.3 Excitation o f T E waves From eq. (19), it is seen that TEm0 wave has the sole electric component Ey given by Ey = Ey o sinAx e~^ (1-25) where Ey o is a constant If the origin of the coordinate system is placed at the center of waveguide’s opening, we have Ey =EyQ cosAx e~y‘ (1-26) For lossless dielectric, y2 = _/'pz, where Pz is a real constant, and eq. (26) becomes Ey = —^-[exp(-yAx - j $ zz) + expfjAx - ypzz)] (1-27) Furthermore, it follows from eq. (21) that, p 2= c o V so eq. (27) can be = > i2 + P z C1'2 8 ) rewritten as £ Ey = —^-[exp(-jP(nrsin0 + zcosG)) + exp(-jP(-xsin8 + zcosQ))] (1-29) where 0 = ±sin_1(-^-) = ±sin_1( - ^ - ) P flP (1-30) Eq. (29) represents two uniform plane waves propagating in the dielectric medium with propagation constant P and direction angles 9 given by eq. (30). Thus, any TE^o wave can be generated by exciting two plane waves propagating in the given directions. These plane waves can be excited by some incident field coming from outside the dielectric medium as in practical microwave heating configurations. Consider a plane wave incident upon the air-dielectric boundary at angle 9,- shown in Fig. 1-5. The plane wave enters the dielectric medium at angle 9 determined by Snell’s law R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited without perm ission. 17 Fig. 1-5 sin0 ; c v1 sine =n = 7 7 = aj7p ( 1 ' 3 1 ) where v \ and v 2 are the velocity of plane wave in air and dielectric, respectively. From eqs. (30) and (31), we get - . A c mn cm 7t sinG; = — ——-— = c o /(3 p a ,, ^ (1 -32) co a To generate a particular TEOTo mode, the incident angle 0,- is thus not determined by the permittivity but only by the dimension of the container and the frequency. Strictly speak ing, Snell’s law for ray optics only applies for an infinitely extended boundary. For a finite boundary as in the actual case of a rectangular container, the refracted field is not only the plane wave with refraction angle 0 but rather a spectrum of plane waves entering the dielectric medium at different angles [14]. For lossy dielectrics, a TE,„o wave can still be expressed as sum of two plane waves. However, these plane waves are non-uniform, that is, the field magnitudes are not uniform but exponentially decaying in one direction along the wave fronts. The refrac tion angle 0 is no longer given by eq. (30) but by a more complicated expression [15]. For most practical lossy dielectrics with relatively small loss tangent, however, eq. (30) R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 18 still gives a good approximation for 0 . Despite these limitations, as we will see in Chapter 2, the excitation method dis cussed above is fairly effective in generating a given mode in a given container although it is in general more difficult to generate a pure high order mode than to do so for a lower one. Finally, it can be shown that any T E ^, mode can be expressed as sum of four plane waves and hence can be generated in the same manner as above. R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. 19 2. WIRE GRID MODELING OF CONDUCTING AND DIELECTRIC BODIES The idea o f using a wire grid model to represent a dielectric body has its origin in the established concept that a continuous perfectly conducting surface can be success fully modeled using a surface wire grid [1]. When illuminated by some electromagnetic source, the scattered field produced by the wire grid approaches the true field produced by the continuous surface if the wire grid spacing is chosen to be small enough. Since no field or current can exist inside a solid perfectly conducting body, it is electromagnetically equivalent to an enclosed perfectly conducting surface and hence can be modeled approximately as a surface wire grid. On the other hand, if the body is dielectric, fields and currents vtithin the body are nonzero which in turn necessitates a volume wire grid model for the body as will be discussed in this chapter. 2.1 Modeling o f continuous media using wire grids Consider the two-dimensional case where a perfectly conducting plate, shown in Fig. 2 -la, is excited by some external source. The induced surface current is in general a vector quantity with x- and y-components, so a wire structure representing the plate must have both x- and y-oriented wire segments, as shown in Fig. 2 -lb. A rectangular patch of the plate is thus represented by a rectangular wire loop of the wire grid. As mentioned before, the choice of wire grid spacing d and wire radius r has led to many discussions in the past [9 ]-[ll], [30]. Intuitively, one can expect that the smaller the wire grid spacing becomes, the better the wire grid should represent the continuous surface. However, in practical modeling, the number of wire segments must be limited and hence there is a maximum allowable wire grid spacing at which the wire grid will still represent the sur face satisfactorily. It has been generally found that d = \I 10, where X is the free space wavelength (at the actual frequency), is adequate to give good results [11], [30]. In many cases, however, a larger wire grid spacing can be tolerated, practically up to X/4 in the R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 20 b) a Fig. 2-1 P-rfecly conducting plate and wire grid model cases where the current distribution is not undergoing quick changes throughout the sur face. On the other hand, no definite answer has been given regarding the choice of the wire radius. Most authors, however, agree that a radius o f r ~ d H k would be the best choice. Ludwig [11] suggests that the rule r —d t2 k must always oe followed and he shows quite cleariy that too large a radius is just as bad as a too small a radius. Consider the solid three dimensional dielectric body shown in Fig. 2-2a. The body can be thought of as being made up o f elementary dielectric cubes, one of which is out lined in the lower left part of the body. When an electromagnetic field is incident upon the body, the induced volume currents require a volumetric wire structure having x-, yand z-oriented segments: such a structure with equal-length segments is shown in Fig. 22b, in which part of the structure outlines tne elementary cube. At mis stage, however, it R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 21 is still unclear how the wire grid should be modified to model the body with a given per mittivity. To answer this question, let us look at the logical build-up o f the model as fol lows. Consider an elementary dielectric cube shown in Fig. 2-3a, which can be represented by the elementary wire cube shown in Fig. 2-3b. Since the body is in general a lossy dielectric, conduction and displacement currents actually flow within the body and the dielectric cube cannot be correctly represented by the wire cube unless it is loaded with lumped impedances to take into account the conduction and displacement currents. The idea of using lumped impedances to represent the loss and phase shift asso ciated with these currents is similar to using lumped impedances to model a transmission b) a) Fig. 2-2 Dielectric body and wire grid model d a) b) Fig. 2-3 Elementary dielectric cube and its wire cube model R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 22 line approximately. If the material has the complex permittivity e = (er - j e r )eo, the effective conduc tivity is defined by a e^ = E r, e 0 CL) (2 -1 ) Furthermore, if the x-, y- and z-components of the volume current are uniform over the respective cross sections, the resistance related to a current component, say Ix, is R= d <5efjA 1 Geffd 1 £ r EgOki (2-2 ) Similarly, the capacitance is given by C = ( e ; - l ) e o — = (er - D M a (2-3) Note that the elementary wire grid is physically located in free space with relative per mittivity equal to unity in comparison with er of the real medium so, the "additional" capacitance must be given as above. The resistance and capacitance must be distributed to the four segments oriented in the x direction, which gives a total load for each wire segment to be 4R in parallel with C /4. In the thin-wire program, it is only possible to load a segment at the ends and in order to avoid asymmetry, the total lumped impedances must be distributed equally at the two ends of each segment as shown in Fig. 2-4. It follows that 4R HH C/4 R. HH C, Fig. 2-4 Total load and equivalent end Lads for an edge segment R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission. Re HH C, 23 (2-4) and (2-5) The above resistance and capacitance are connected in parallel because the conduction and displacement currents add in Maxwell’s equations. Thus a lumped impedance is gen erated and this impedance is inserted at each end of each segment in the elementary wire cube. The dielectric body shown in Fig. 2-2a consists of five elementary dielectric cubes and it can be modeled by joining five elementary wire cubes together as shown in Fig. 25. Consequently, groups o f two, three or four segments o f adjacent elementary wire grids are merged into single segments which we will call surface -, comer - and volume seg ments, respectively, corresponding to their positions in the resultant model. Those seg ments that are not merged with any other segments are termed edge segments. The lumped loads for the first three kinds of segment are the parallel combinations of the edge impedances given in eqs. (4) and (5). Of course, one can only tell the differences among these kinds of segment by looking at their respective lumped loads. Without the lumped loads, all segments will be physically identical. In terms of edge impedances, the surface -, comer - and volume impedances are given by ( 2 -6) a dielectric body is thus completely represented by a wire grid model loaded with lumped impedances, where the wire segments allow calculation of scattered field from R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 24 / \ Surface segment w ith/?,, C, Edge segment with/?,, C, Com er segment w ith Rc, Cc Volume segment w ith/?,, Cy Fig. 2-5 Building up the wire grid model of the dielectric body shown in Fig. 2-2a by joining elementary wire cubes. the body and the lumped impedances allow calculation of loss and phase shift within the body. We will concentrate on modeling bodies using cubic elementary wire grids in this thesis. Non-cubic but still rectangular building blocks can be modeled easily following the same principles. Further study of non-rectangular building blocks is required before bodies of arbitrary shape can be modeled. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission . 25 (0.5,0.5, 1.0) — source dipole 6 / ( 0 .01 , 0.01 , 0 .01 ) / ■ • / / center / / / ?. = 0.12 m at / = 2.45 GHz (0.00, 0.00, 0.00) Fig. 2 -6 -10 -15 -20 tu S. f=2.45 Ghz -25 iu d=10 mm -30 -35 ^0 -45 0.4 0.6 0.8 2.4 1.2 r[mm] Fig. 2-7 Normalized field magnitude vs. wire radius R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission. 26 Before going over to computation with various models, let us try to confirm the r -dl2.il rule as the best choice of wire radius in volume wire grid modeling. It should be noted that an unloaded elementary wire cube represents a perfectly conducting block, within which the field should be equal to zero. Consider the configuration shown in Fig. 2-6 where the wire cube is illuminated by the sinusoidal field from a source dipole placed some distance away. The field magnitude at the center of the wire cube is computed and plotted against the wire radius as shown in Fig. 2-7. Note that the field has been normal ized to the incident field, that is, the field that would exist without the wire cube, to give some idea of the effective shielding of the grid. The normalized field takes a minimum at r = 0.0012 m, which gives dir = 8.33 which is somewhat larger than the expected ratio of 2k. Aversa [3] conducted a similar experiment with a larger hollow box and he found the best spacing-to-radius ratio to be about 7. Since the radius step in his experiment was big, the minimum could as well be somewhat larger or smaller. However, doubling or halving the "best" radius definitely gives worse results. As a conclusion about the choice of wire radius, the rule d =2kr should be followed but only approximately. In general, the ratio d/r = 2 k will be used through the modeling work contained in this thesis. 2.2. Dielectric heating in a small, deep rectangular container Having established the fundamentals for wire grid modeling of materials, we wish to validate the approach by modeling and computing the field distribution within some rectangular containers filled with lossy dielectrics. Besides having the practical value of objects that could be used in microwave heating, rectangular containers seem to offer a particularly useful way to validate the method through comparison of numerical results with analytical solutions known from the theory of rectangular waveguides. Consider the rectangular container filled with lossy dielectric shown in Fig. 2-8. The container must be chosen so that only the lowest modes can propagate while higher modes are suppressed. Ideally, we wish to have only one mode within the container so R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. air-dielectric boundary b/2. perfectly conducting container z Fig. 2-8 Geometry of container filled with dielectric that the analytical solution is exactly known and can easily be compared with the result deduced from the wire grid modeling approach. It is generally known from the theory of waveguides that waveguides with sufficiently small cross section can support only the lowest TEjo mode, while being cut off for the higher modes. This can be seen from eq. (1-27), where the propagation constant yz is essentially a real constant for high order modes, which means these modes attenuate rather than propagate along the z direction. From this point o f view, it is desired to choose a and b as small as possible; a should be large enough to support only the T E 10 mode and b can be chosen arbitrarily small. Mathematically, we have v - I, V a s?— TTrnr?----T T -’— "7 T b R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 28 U b < a, the cutoff condition for all modes except TEio requires ( - ) 2 - cd2|1£^.£o < 0 < (— )2 - c o V ; e o a a (2-8) or k 2k — — — <a< — — — wV(i£r eo “ VliCrEo (2 -9 ) For butter with er = 3.5—j 0.9 at / = 2.45 GHz, a can thus be chosen in the interval 0.0327 m < a < 0.0654 m. We choose, however, a = 0.07 m in order to have more data points for the evaluation of the numerical result. With this value o f a, the TE 30 mode will be cut off but the TE 20 mode will not. By using symmetrical excitation, however, all odd TEmo modes, that is TE2o, TE 40, ... will be suppressed and the dominance of the TE 10 mode is ensured. The depth I of the container must be large enough to show the field variation along the z direction. For the TE 10 mode, the propagation constant is deter mined from eq. (7) with m — 1, n = 0 to be yz = 13.79 + j 86.05 m-1, which gives the wavelength of the TE 10 wave Xz — 0.073 m. We choose / = 0.05 m. The wavelength o f a plane wave in an infinite dielectric medium is given by X = y , (5 = l m ( y } = I m { ^ J I ^ } (2 - 10 ) which gives X = 0.065 m. The grid spacing can be chosen to be d = 0.005 m, which is about X/13. Now, the last dimension b can be chosen as small as 0.02 m. In fact, a smaller value of b could be chosen but this would not allow us to see the field variation along the y direction, which is necessary to determine the existence of certain modes within the container. With the chosen values of a, b and /, the size of the wire grid model for the container becomes (14x4x1 OK As the origin of the new coordinate system is placed at the center of the container’s opening, the container and the filling dielectric material are fully symmetric about the x and y axes, so that electric and magnetic sym metries can be utilized to reduce the number of wire segments and dipole current modes. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission. 29 The fundamental theory of these symmetries is treated in Appendix A. Hence, only a quarter o f the full wire grid model needs to be set up. Fig. 2-9 shows the quarter wire grid models for the container and the filling dielectric separately. Note that the segments at the edges, on the surfaces and inside the dielectric model are loaded with edge-, surfaceand volume impedances respectively and are determined from eqs. (4)-(6): Re = 3260 Q. Ce = 0.55 pF Rs = 1630 Q Cs = 1.10 pF Rv = 815 f t Cv = 2.20 pF The effect of enclosing the dielectric with a perfectly conducting container is that the lumped impedances on the segments shared by both models are shorted. This is a direct consequence o f the fact that wire segments are themselves perfectly conducting and as building elements for the container, they are unloaded. The complete dielectric-in container wire grid model is shown in Fig. 2-10. Also shown in this figure is the source monopole oriented in the y direction. Using ME-symmetry in the XY-plane, the source configuration becomes two dipoles placed symmetrically about the X = 0 plane at a dis tance R = lm from the center of the container’s opening. This distance is large enough so that the radiation field from the dipoles approximates two y-polarized plane waves incident at angle 0,- to the z axis. To give a proper excitation for the T E ^ mode, the incident angle 0/ is chosen according to eq. (1-38): „ . _i,cm it, . _i, 3T08T n_ , c , o -------- ) = sin ------------ 3------- ) - 6 1 co a 2-2.45T0 -0.07 0 , = sin The wire grid model of a dielectric in a container together with the excitation source form a complete wire configuration as input for the thin- wire program. A fter the induced currents on each segment have been computed, the near field at a particular location can be determined as the sum of the radiating fields from all the segments. In general, the field values are most reasonably sampled at the centers of the wire cubes since these are R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 30 // t 4 V/ / f/ A/ f,/ V/ V* V/ Vt ■ dielectric - T /T 7 7 7 7 L a L 1 7 i i l J z. 7 i \ I 2 ? 2 i 2 2 I 7 2 7_7 i I 2 2 7 \ 2I I I 27 7 2I I 7 17 7 t I i. 2 I 2 7 7 I I 7 77 AA 2 t t 1I I l ? 7 L AI 2 1 I 2I 2 l ?7 LA 2 I I AC A 1 1 1 J 2I J 2 7 7 r t t * / / / / / / / / I I II I I II I I I I I I? I I 5 .A 1 E EEEE /.J container d=0.005 m a=0.07 m b=0.02 m 1=0.05 m Fig. 2-9 Quarter wire grid model for container and dielectric source monopole R 6, M a/2/ f - - 0 ; = 61 ° IL ACr ' a <- AL i l l t j l V- I I I I I I Q l in H l V J 11 11 ■H R=1 m IV V Fig. 2-10 Complete quarter wire grid model for dielectric in container R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 31 the building blocks of the model. Fig. 2-11 shows the variation of the \Ey | components along the x axis for y = 0.0025 m and z = 0.0025m; the field values are normalized to the value at y = 0 in order to be compared easily with the cosine distribution shown in the solid curve. The numerical results follow the cosine distribution closely, indicating that no higher order modes of m > 2 are present in the actual field. It follows, therefore, that only the TE10, TEH and TM U modes are contributing to the total field. Fig. 2-12 shows | Ey | and | Ez | as functions of x for different values of y and z = 0.0025 m. The Ex com ponent is negligible and therefore not plotted. Clearly, the Ey is the dominant component, indicating that only a small portion of the TMu mode is present. Furthermore, the small variation o f Ey along the y direction shows that the T E n mode is also insignificant. Together, these modes represent no more than 1/10 of the total Ey component. Fig. 2-13 shows the variation of normalized Ey component along the z direction. It appears that the numerical result displays a stronger attenuation than the theoretical standing wave pattern for the TEio mode shown in the dashed curve, although the deviation is not large. A mix ing o f TEio and T E n modes by the ratio 10:1 and phase difference <{>= 7t produces a pat tern matching much better with the numerical result. This is no surprise since the T E n mode is cut off and hence attenuates more than the TEio mode. A contribution of the T E U mode to the total field is therefore likely to intensify the attenuation. The numerical result obtained is in fact very good since the propagation and the attenuation presented in the standing wave pattern is very sensitive to any change of the complex pennittivity e. Tests have been done with other dielectric materials and the numerical results agree quite closely with the analytical ones. Power dissipation in each wire cube can be determined by simply summing the ohmic losses on the lumped resistances belonging to the ware cube P = Z T ( «i l / i I +*«2/?2 ) (2-11) i=l 1 where —R[\ l}\ and —rt /2 are the ohmic losses at the two ends of i ’th segment of the R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. 32 0.9 lEyl/lEyOl 0.7 0.6 0.5 : cos(x*pi/0.07) 0.4 0.3 0.2 (x): numerical result y =0.0025m , z=0.0025m 0.005 0.01 0.02 0.015 0.025 0.03 0.035 X [m] Fig. 2-11 Normalized field distribution Ey at the opening of rectangular container filled with lossy dielectric 50 45 Field Magnitude (V/m) 40 35 30 25 20 15 10 (x): lEyl for y=O.025m. z=0.025m (o): lEyl for y=0.075m, z=0.025m (*): lEzl for y=0.025m. z=0.025m (+): lEzl for y=0.075m. z=0.025m 5 0 0.005 0.01 0.02 0.015 0.025 0.03 0.035 X [ml Fig. 2-12 Distributions o f some field components in container R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission . 33 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 - : Analytical solution : TEIO : Analytical solution : TEIO + 0.1*TE11 (o): Numerical result for v=0.0025m , z=9.0025m Z 1m Fig. 2-13 Standing wave pattern in rectangular container filled with lossy dielectric er = 3.5 —j 0.9 y ^ .0 0 2 5 m , r=0.0025ra 0.9 0.8 (x): numerical result 0.7 : (cos(x*pt/0.07)) 0.6 o 0.5 0.4 03 0.2 0.1 0.005 0.01 0.02 0.015 0.025 0.03 0.035 X [ml Fig. 2-14 Power distribution at the opening of the container I I \ R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 34 block, and c, = 1, 1/2, 1/3 and 1/4 for edge-, comer-, surface-, and volume segments respectively. This factor is to take into account the ohmic losses shared by adjacent wire cubes with common segments and lumped resistances. Fig. 2-14 shows the distribution of normalized power dissipation along the x axis. The numerical result agrees well with the theoretical cosine-squared distribution although slight deviations are observed for small values of PIP o- The most probable explanation is that the contribution of the Ez component to the total power dissipation becomes apparent when the dominant component Ey is small. There is likewise good agreement between numerical ar.d analytical results for the power distribution along the z axis but the result is not shown. 2.3. Dielectric heating in o th e r rectangular containers The container considered in the previous section is only a special case, in which the dimensions have been chosen together with a special source configuration to support only the lowest TE mode. In practice, containers can have various shapes and sizes, and are usually exposed to a very intricate microwave field such as the one in a commercial microwave oven, which can have many polarizations and magnitude distributions. The power dissipation in the dielectric can therefore have a very complicated pattern due to the complex field distribution within the container. Consider for example the container and source dipoles shown in Fig. 2-15. The dimensions of the container have been chosen to support several modes beyond the lowest one. The power density distributions in the XY plane for different z values are shown in Fig. 2-16. It is noted that the plots are of different scales to demonstrate the shape variation of the power absorption along the z axis. The power distribution at z = 0.25 cm is seen to have two main lobes at x = 0 and four smaller lobes at the comers. It appears the shape is formed by a mixing of T E ^ , TE/TMU , TE/TM12) TE/TM 31 and TE/TM 31 modes with TE/TM n being the dominant R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission. 35 2cm source dipole 10cm E -s- lm lm 7cm .5cm Container Dielectric fU ^ rv J J Fig. 2-15 Heating of dielectric (er=3.5-/0.9) in (20xl4x3)d container with d = 0.5 cm being the grid spacing for the model Fig. 2-16 Power density distributions at different depths inside the dielectric in container shown in Fig. 2-15 R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 36 components. As different modes have different propagation and attenuation constants, the shape of the power distribution changes along the z axis. Symmetries are obtained about the x and y axes because only one source configuration ME is used in this case. Non-symmetric power distributions may result by combining different source configurations (see App. B). One of the motivations behind microwave heating research is to obtain heating uni formity in food materials. If one could generate any mode with desired magnitude and phase in a given container, it would be possible to synthesize a uniform heating pattern by superimposing different modes with properly chosen magnitudes and phases. The difficulty is, o f course, how to generate these modes. Even for a container large enough to support a certain high order mode, it is difficult to control the mode excitation because of the diffraction over the edges of the container. Consider the container shown in Fig. 2-17. The model and source configuration is similar to that shown in Fig. 2-10. The source dipoles are at distance i? = l m from the center of the container’s opening and several incident angles 0,- are chosen in the interval [0°,90°]. The heating patterns along the x axis at y = 0.00375 m and z = 0.00375 m for different incidence angles are shown in Fig. 2-18. It should be noted that 0/ = 17° and 0,- = 61° are incidence angles for TEio and TE 30 modes, respectively, according to eq. (1-38). The heating patterns for these angles exhibit approximately single mode field distributions as expected. At other angles, the heating patterns show a mixing of T E ^ and TE 30 modes with different amplitude ratios and relative phases. It is interesting to note that for 0X - = 0°, the contribution of the TE 30 mode is quite significant while it completely vanishes at 0,- = 17°. For 0; = 90°, the power absorption in the dielectric arises totally from diffraction of the incident field over the edges of the container. Further study of field transmission in containers made of con ducting and dielectric materials seems necessary to optimize the heat distribution. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 37 M 2cm source dipole R=lm ' 21cm E 3cm 2.25cm dielectric Fig. 2-17 Heating of dielectric (er=3.5- j 0.9) in (20x4x3)d container with d = 0.75 cm being the grid spacing for the model 0.02 0.018 0.016 0.014 m 0.012 E •Si a. 0.01 0.008 0.006 0.004 0.002 -0.15 - 0.1 0.05 -0.05 0.15 X [m| Fig. 2-18 Power density distribution inside the container shown in Fig. 2-17 for different incidence angles. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 38 3. ONE-DIM ENSIONAL M ICROW AVE HEATING - TH ERM A L RUNAWAY Non-uniform distribution of heat in dielectric materials when processed by microwave power may be caused by the limited penetrability of the electromagnetic wave as well as by standing wave or other effects. This non-uniformity is sometimes intensified by the inhomogeneity of the material where the above effects are reinforced. Due to the temperature dependence o f the relative permittivity, the heat generated within the material can change the dielectric structure in such a manner as to concentrate the power dissipation more and more into the hottest regions, which leads to thermal runa way with damaging effects for the material. Besides the problem of determining the field distribution in an inhomogeneous material, the mutual dependence of the power dissipa tion and temperature through the dielectric properties obviously contributes to the com plexity of the microwave heating analysis. It seems therefore natural first to investigate the problem of one-dimensional microwave heating, which is probably an over simplified case but nevertheless will give us some understanding and insight into the phenomenon of thermal runaway. 3.1. W ave propagation in inhomogeneous media Consider a uniform plane wave incident upon the surface o f a dielectric half-space where z > 0 as shown in Fig. 3-1. The relative complex permittivity o f the dielectric medium is assumed to be independent of x and y but is a function of z, i.e. z(z) = (er( z ) - j z ' r (z))£o (3-1) Also, assume the only component o f the electric field is E = Ey . Since there is no field variation in the x-y plane ( -^ - = - ~ = 0 ), the wave equation given by eq. (1-13) is ox dy reduced to R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 39 E(2) eo E: propagation z=0 Fig. 3-1 Geometry for one dimensional microwave hearing d 2£ 2 - y 2(z ) £ ' = 0 y2(z) = -co2 ^ 0 e(z) (3-2) This equation for wave propagation in inhomogeneous media can be solved by various approximation methods. The conventional WKB method is most convenient for finding the solution for slowly varying media. Let E(z) = A(z)e~i m (3-3) where A and 4> can be regarded as the amplitude and phase, respectively. Substitution into the wave equation gives A" - 2y<t>'A' - ;<(>"A + [—"y2 (z) - ((j)')2]A = 0 (3-4) where the prime indicates a derivative with respect to z. For slowly varying media, we make the approximation a" = 0 Then A e ^ isan approximate solution of eq. (4) if <j> = ± 7 t(z ) or 0 = ± yJ t z ) d z and R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. (3-5) 40 or A= const (3-6) ■ fT W f then, the general solution is a exp - J V E(z) = (z ) dz + b exp + J y ( z ) dz 'F J W ) (3-7) o where a and b are arbitrary constants. Let y(z) = a(z) + y‘P(z), the integrals in eq. (3-7) representing the propagation and attenuation of two waves travelling in the positive and negative z directions. Since the medium is infinite for z > 0, only the first wave exists. Furthermore, additional amplitude and phase variations of this wave are given by the fac tor outside the brackets. The field amplitude is given by 2 \E(z)\ = V |i(z )-y a (z ) exp - J a (z)dz (3-8) d$ It is interesting to note that for a low-loss dielectric, a c (J and if - j - < 0, we could have a growing wave even in a lossy medium, at least at some short distance from the surface where the decaying characteristic of the factor in the denominator is dominant over that of the exponential function. On the other hand, if dz > 0, then we will have another attenuation factor beyond the one caused by the loss of the dielectric. This attenuation can be regarded as being a result of the internal reflection at every point along the z direction caused by the growing characteristic of the permittivity. The generalized WKB method is well described by Wait [16], who gives a more exact solution for the wave equation. Anyway, even this modified method must be based on the assumption of slowly varying media. For rapidly varying media, other methods are available [17], but they are less compact and not necessarily numerically efficient A numerical approach based on the finite difference method will be utilized in the following sections to solve practical problems of one-dimensional microwave heating. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 41 3.2. Numerical method for coupled wave-thermal equation Microwave heating is usually a dynamic process, in which most of the parameters are functions of space and time. The field distribution in the dielectric medium has a direct impact on the temperature distribution through the power dissipation given in eq. (1-16). On the other hand, the field distribution depends on the dielectric profile, which is a function of the temperature distribution. The time dependence of the electric field mag nitude and dielectric properties comes into the picture as the temperature itself is time dependent. The transient electric field need not be considered since it practically will not have any impact on the slow thermal process. The one-dimensional thermal equation is given by dT d2T p 3or7 = a'T9 zTz + TT" pcp (3-9) where T = T ( z , t ) : temperature [°C] p = p (z,r) : power density [W/m3] ~ cq = ----- : thermal diffusivity [m /r] PCp Xc : thermal conductivity [W/(m°C)] p : mass density [kg/m3] cp : specific heat [J/(kg°C)] The first term on the right hand side of eq. (9) represents the rate of change of tempera ture by thermalconduction and the second one represents the heatgeneration. Normally, thermal conduction can be ignored due to the quick microwave heating process. Since eq. (9) presents no major problem when solving numerically, thermal conduction will be included for the sake o f generality. Assume the temperature dependence of the relative dielectric constant and loss fac tor is given by R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 42 Er( z , t ) = f ( T ( z , t ) ) and er = g ( J (z,r)) (3-10) Combining eqs. (2),(9),(10) and (1-16), we get the coupled wave-thermal equations dz + co2 Moeo[f(T) - jg (D J£ = 0 (3-11) (3-12) 3.2.1 Boundary condition considerations Eqs. (11) and (12) are second-order differential equations, which can be solved with tv/o specified boundary conditions for each equation. First, consider the wave equa tion for the half space problem shown in Fig. 3-1. Theoretically, the wave equation can be solved by the Runge-Kutta method [18] when the boundary conditions £ (0 ,r) and E (0,r) are given. For a homogeneous medium (for which the exact solution is known), E (0 ,0 and £ (0,0 can be found if £,• is given, as £(0,/) = r£f (3-13) where T is the reflection coefficient at the surface. £ (0 ,r) = - y £ ( 0 ,r) (3-14) The second boundary condition is valid as only a decaying wave is assumed. Of course, there is no reason to solve the equation numerically when its analytical solution is already known, but this discussion will outline the problems associated with the wave equation. With the given boundary conditions above, the wave equation cannot be solved in practice because of numerical instability. Recall that the complete solution for the wave equation is a linear combination of decaying and growing waves. With the boun dary condition given in eq. (14), only a decaying wave should be selected but the inevit able round-off errors in the numerical solving procedure will excite the growing wave. Eventually, this will become the dominant component causing an erroneous solution with wild oscillations. Alternatively, one can use the finite difference method, implicit scheme R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 43 [19], to suppress the instability but the solution is of poor accuracy. At this point, it seems that any numerical method using the above boundary conditions cannot yield a satisfactory solution. On the other hand, as the actual medium is lossy dielectric, the field at infinity vanishes; therefore, practically, the field is negligible at some distance z = d, which can be estimated as follows (3-15) which gives (3-16) where 8 is some small constant and a min= m m [a(z)] 0 < z< °° (3-17) Intuitively, the attenuation in any medium with a(z) > a , ^ for all z must be larger than that in a homogeneous medium with a(z) = for all z. Thus, eq. (16) gives the dis tance at which a vanishing field is ensured. With the boundary conditions E (0,r) and E( d, t ) = 0, the wave equation (11) can be solved using the finite difference method and the solution is guaranteed to be stable with good accuracy. The thermal equation (12) also can be solved using the finite difference method with the initial conditions (3-18) 7(z, 0)=7~o(z) and the boundary conditions for the thermal isolation case r'(0 ,r) = 0 and T{d, t) = 0 (3-19) or, for the case of convection loss at the surface z = 0 T (0,t) = -^-[Tan - 7 (0 ,0 ] and T(d,t) = 0 (3-20) where h is the convection heat coefficient [W/(m 2 oC)]. In both cases, the boundary R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 44 condition T (d,t) = 0 follows from the negligible field and hence power dissipation at z = d, which means a vanishing temperature differential at that location. 3.2.2 Numerical approach The space-time domain can be discretized into a finite difference grid as shown in Fig. 3-2. Suppose the solution domain is restricted to (M x N ) nodes. Then Xnun= X [ ( m - \ ) h , ( n - l ) k ] \< m < M ,l< n< N (3-21) where X is any quantity varying with z and t. Wave Equation. The finite difference formula for the second order derivative of E at node (m ,n ), n constant, takes the form d~E Em+\,n "h Em —l.n ~ 1<m <M — —= ---------------dz2 h2 (3-22) With substitution of eq. (22) into eq. (11), the wave equation takes the finite ditference form at time step n (■ •i ■>■— ■ »— (> ■ <i ■ — <> — — <3 f < (m.n) ( k 1r ► • • • . m= 1 x h m==M Fig. 3-2 Finite difference grid of space time domain R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission . 45 E m -l,n ^ m ,n + E-m+l.r. ~ 0 ( . 2 + h \ < m < M (3-23) where T ^ n = “ 2 ^ e o [ /( 7 ’m,n)-;'g(7'/n,n)] (3-24) Assuming £ l n and E m „ are given as boundary conditions, wehave from (23) a system of (M - 2) equations,which can be solved easily. Usually, only theincidentfield £,• is given and E 1>n = E (0) must be determined from eq. (13) with t(°) = — n(°)~ ■ Ho+ 11(0 ) (3-25) where TJo = yPoTSo" 2nd r](0) is unknown. It should be noted that, by definition n(z) = E(z) H(z) E(z) 1 E ' {z) J'mo For homogeneous media, E / E = y and hence rj = (3-26) qj/e. For inhomogeneous media, this is no longer valid since the solution is still unknown. Instead, we assume E(z) - C E(z) (3-27) where C is an arbitrary constant. We start to solve the boundary value problem r 32£ 3z 2 te =0 £(0) = C E (0) and E(d) = 0 (3-28) Once £(z) has been found, its derivative can be determined by numerical differentiation [20] E(Zs) = I ^ - E ( z u) + R [ E ] u=0 h where R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. (3-29) 46 , if S Lu V_1 = U i v v*u (3-30) -1 I I , if 5 1 1 (Z|» —Zv) u v=l V= 1 v*u v^s.u For five-point interpolation (/ = 5) we have E (0) = — [-25£(0) + 48£(A) - 36E{2h) + 16£(3/r) - 3£(4/t)] + 4-A4£ (5)(3-31) 12/i It readily follows from eqs. (26) and (27) that tKO)=7O)Po-I^0) (3-32) £ (0) E (0) is then determined from eqs. (13) and (25) and hence E (z) can be found by E ( z) = ! ^ C = - P ^ £ ( z) E( 0) (3-33) Thermal Equation. Equation (12) involves both space and time variables for which a widely used differencing scheme of Ciank and Nicolson can be found in [21] T m,n + 1 ~ T 'm ,n Ctf r k l 7--------- = ~ [ T " >• . 1 m,n + F m,„+l] + —[Gm,n + Gr^n+i] 2 ^ ^ 0 < m < M (3-34) where r m + l,n + T m - l.n “ 2 F m_„ nun _ Gm n coeo pcp 2 I^-m ,n I Eqn. (34) can be rearranged into the form R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. (3-35) (3-36) 47 - j - T m. hn+l + (1 + r)Tmn+i - y 7 m+1.n+1 = Dn 0 <>m<M (3‘37) where r = -Kr and h2 = [~2^'m+l,n ~ (1 —r )Tm,n + + “ [^m,n + Gm,n + 1] (3-38) The boundary conditions given in eqs. (19) at the time step (n+1) take the finite differ ence forms T'i.n+i -T i,n + \ and TM'H+i = r^ /_ l n+1 (3-39) and for eq. (2 0 ) (1 + “ ) r 1.„+1 = r Zn+ 1 + - ^ r am and r M.fl+1 = r lW_1.n+1 (3-40) From eqs. (37) and (39) or (40), we have again a system of equations with ( M - 2) unknowns. These are the temperature at time step n+ l expressed explicitly by constants and quantities known from time step n as given in Dn, with one exception of GOTin+1. The coupled wave-thermal equation can be solved iteratively as follows: 1/ At time step n with known temperature distribution Tmn, solve eq. (23) for field distribution E ^ ni 0 <m <M. 2/ Find G ^ „ 1 <m <M from eq. (36) 3/ Assuming Gm,n+1 = Gm,„, solve eq. (37) for T„^n+X 4/ Repeat 1/ and 2/ for time step n+1; solve for Tm^n+\ with the new value of . This procedure is repeated until Tm<n+i converges. Usually, the conver gence is robust and one iteration is enough. 5/ Go to 1/ for next time step. It is noted that both system equations (23) and (37) are of tri-diagonal form. Using Gauss elimination, the solutions can be expressed in implicit formulas and then found by back substitution [2 2 ]. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 48 3.3. T herm al runaw ay in one-dim ensional microwave heating Thermal runaway can be defined as a process -a which the difference between the n.tes of temperature rise at two particular points grows with time, or mathematically d \ T ( z a) - f ( z b)\ (3-41) 0 dt where the dot indicates a derivative with respect to time. This definition is consistent with the common understanding of thermal runaway as a process in which the power dis sipation becomes more and more localized it some particular points. The relation is obvi ous from eq. (9) by neglect'ng the conduction term T(z) = E llL =:G( z ) Pcp (3-42) which can reasonably be done as thermal conduction is usually small in microwave heat ing. G (z) will be referred to as heat source density. Substitute eq. (42) into (41), we get d \ G ( z a) - G ( z b)\ dt (3-43) >0 Thermal runaway can thus be illustrated graphically as in Fig. 3-3 which sketches heat source density distributions at two different times t\ and r?, where t ] < t 2 a) b) a O z z Fig. 3-3 Postulated qualitative distributions of heat source density for cases of a) thermal runaway and b) non-thermal runaway R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 49 3.3.1 Examnle 1: Cooked beef juice The dielectric properties of cooked beef juice at / =2.45 GHz [6 ] are approximately given as functions of temperature T er(T) = 7 0 - j T e” (T) = 23 in the temperature interval [5°C , 65°C]. Assuming the thermal properties of cooked beef juice are the same as those of water, we have cp = 4180 J/(kg°C), X, W/(m °C), p = 1000 kg/m 3 and h = 6 W/(m2 CC) for free convection with air which causes a negli gible surface loss as seen bv eq. (20). Assume a sinusoidal field of magnitude t i = 20000 Vh n incident upon the surface of the half space dielectric medium as shown in Fig. 3-1. The distance at which the field magnitude E(d) < 10- 3£ (0 ) is determined from eq. (16) is ^ _ InS ____________________ Int.iO-3)___________________ = 0 1 m c min " Re(-(2rt 2.45T0 9 )24 k - 10-78.85H 0 -12(70 - j 23)} 'm Suppose the initial temperature distribution is T (?. )) = 10 °C f o r C < z < J . With the numerical method developed in the previous section, a simulation o f 1-D microwave heating process is run ovei a time period of rmax = 8 .s. Tne length and time ster are h = 0.0005 m and k = 1 s, respectively. The temperature and power density distributions at different elapsed times are shown in Fig. 3-4a and b respectively. As the loss factor is constant, the field distribution is directly proportional to the square root of the power density and is not shown. The results obtained show a weak thermal tunaway taking place at the surface of the dielectric medium: the heat source density at the surface increases from 6.6 to 7.8°C/s while it remains unchanged for z >0.01 m. The explanation to this observation, despite the constant loss factor er = 23, can be sought in the decrement of er at the surface where the heat source density is highest and hence the temperature rises fastest. First of all, the R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. so b) a) 65- 6- G 4°C/s! ■r 60- 2- 0.01 0.02 0.03 7. 0.04 0.05 [m] 0 0.01 0.02 0.03 z [m] 0.05 c) : t = 0s : t = Is 40-Ar T [°C] : t = 8s 20- -e- -0 0.01 0.02 0.03 0.04 0.05 z [m] Fig. 2-* Distribution o f a) heat source density b) dielectric constant and c) temperature in cooked beef juice reflection coefficient T drops for decreasing £r , i.e. the transmitted field at the surface increases. This increasing field, however, decays quickly away from the surface because the medium, after the initial heating, becomes inhomogeneous as shown in Fig. 3-4b. Clearly, the increasing attenuation must be due to the growing dielectric constant profile in the direction of propagation, which in turn causes internal reflection within the medium as has been explained in Section 3.1. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 51 3.3.2 Example 2: Nylon 66 Nylon 66 is a typical example of a material that manifests a strong tendency for thermal runaway [23]. Its dielectric properties are shown in Fig. 3-5. An investigation of thermal runaway for Nylon 66 can be carried out as before assuming £,- = 120000 V/m, / = 2.45 GHz, cp = 4180 J/(kg°C), = 10 W/(m °C), p = 1000 kg/m 3 and negligible con vection loss at the surface of the dielectric half space. Because of the low loss factor, d = 8 m is chosen to ensure a vanishing field at that distance.The simulation time rmax is chosen to be 8 s and the steps h =0.005 m, k — Is. The numerical results are shown in Fig. 3-6. Thermal runaway is obvious by the strong shift of G (z) toward the surface. This is mainly caused by the significant increase of e r at the surface, which is brought about as a dsr « result of the positive slope, ------ , of the er vs. temperature response. Since G is propordT tional to er , tius is the positive feedback that typically causes tnermal runaway in the 3 3 3 100 UO 1G 0 160 T C O 0 5 04 0 3 0 2 01 0 T CC) Fig. 3-5 Dielectric properties of Nylon 66 R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 52 a) b) 80 - | 156010- 20- ■0 4 3 7. •0 1 [m] 2 3 4 3 4 z [m] c) d) .15 2 .9 9 -, t r 2.98 - i i 2 .9 7 - 0.05 2.96 0 2 4 3 -0 i 7 [m] 2 z [m] c) — B : / = Os ■x- : t = 4s 100- -A- : t = 8s 80 T [°C] Fig. 3-6 Distribution of a) field magnitude b) heat source density c) dielectric constant d) loss factor and e) temperature in Nylon 66 60 40 - 0 1 2 3 4 z[m ] R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 53 materials that possess this property. As a secondary effect, the rising er leads to increas ing power absorption at the surface and the field attenuation increases causing a drop of G for z > 0.3 m. It is noted that the transmitted field at the surface appears unchanged despite the rapid rise of Er ; this is because er is small in comparision to e r , which almost does not change in the actual temperature interval [40°C, 80°C] ( see Fig. 3-5 ). There is, in fact, a factor that counteracts the runaway effect in this case; it is the decaying profile of both £r (z) and e r (z), which causes the field to grow in the z direction. However, the effect is too weak to balance the strong field attenuation caused by power absorption at the surface. 3.3.3 Example 3: Thawing o f frozen water W hen defrosted, water and water-containing substances like most food materials undergo a phase change near 0°C, at v/hich significant changes o f dielectric as well as thermal properties are observed. As mentioned in Section 1.1, the relatively small er and e r of ice make the penetration depth o f electric field much larger than that in water. In a thawing process, when the ice at the surface of the material melts, the layer of water will absorb most of the microwave energy, preventing further penetration into the material. Of special importance in the thawing process is the significant latent heat of water at 0°C; that is, the he^'. required for the ice-to-water phase change. In order to incorporate this factor into the simulation, it seems most natural to reformulate the thermal equation (9) in terms of enthalphy and temperature 1 dH d2T p — TT + _ cp dt 1 dz2 pcp (3-44) where H is given by \c pT t H = y CpT + La H : enthalphy [J/kg], <0 T> 0 Lc\ latent heat [J/kg] R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. (3' 45) 54 The thermal equation given by eqs. (44) and (45) will be solved by the finite differ ence method as before. For water, La = 3.33T0 5 [J/kg] and unlike the cases studied before, Xt , p and cp can no longer be considered as constants. The temperature dependences of these quantities [5] are here approximated by piecewise linear 0 W t % functions as shown in Fig. 3-7 along with those for s r and £r . In this simulation we choose Ei = 20000 V/m, / = 2.45 GHz, r max = 8 s, h = 0.002 m, k - I s . The numerical results are shown in Fig. 3-8. Looking at the field distributions, we see a smooth curve that dies out at around z = 0.1 m for t - 0. At this starting moment, the medium is homogeneous because the temperature everywhere is —8 °C, so the field distribution is exponential as expected. Because of the high power density, the temperature at the surface rises fastest and so do er and e r . For t = 1 s, the surface temperature is about —1°C and er, e'r have risen so much that the transmitted field has been reduced significantly. On the other hand, e'r and e r' still remain low for z > 0.02m, which gives rise to a local maximum of the field there marking the transition from a medium with high permittivity to another with low permittivity. For r = 2s, the temperature at the sur face and up to a depth of 0.01m has reached 0°C indicating a phase change is tak ing place. The medium in this "mushy" region is a mix of ice and liquid water hav ing the same temperature because all the dissipated power in the region goes to raising the enthalphy. The widening of the mushy region continues until t = 16s, when the surface enthalphy becomes high enough to overcome the latent heat needed to complete the phase change from ice to liquid water. After this moment, the temperature rises quickly in the newly formed layer of water at the surface and the mushy region moves deeper into the medium. At t = 24s, the temperature at the surface is over 30°C while it is still 0°C just 0.005m below the surface. It should be noted that the rate of temperature rise dT/dt is no longer approximately propor tional to G in the mushy region due the latent heat, so G cannot be used as an R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 55 b) a) 60 • 10 d) c) C 5, 3000 4< at 06 2300 0.4 •10 • 10 e) I os a? 012 0.1 •10 40 Fig. 3-7 Variation of a) relative permittivity u) iuss factor c) tiierniai cunducdvity u) specific heal and c) mass density of water on temperature around 0°C R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission . 70 56 b) a) 10 10- G [°C/s] 5 5- 0 0 0.02 0.04 0.06 0.08 1 ■0 0.02 0.04 0.06 z [m] z [m] c) d) 0.08 1 80 60 20- 40 10- 20 ■0 ■0 0.02 0.04 0.06 0.08 1 -0 0.02 0.04 z[m ] 0.06 0.08 z[m ] — B : t = Os c) 1- : r = Is 30- \ — A- : t = 2s — x- : t = 16s 20- — e- : t = 24s Fig. 3-8 Distribution of a) field magnitude b) heat source density c) dielectric constant d) loss factor and e) temperature in water -10 -0 0.02 0.04 0.06 0.08 z [m] R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 57 indicator for the thermal runaway effect as before. However, one can see that the difference between the rates of temperature rise at the surface and at some point in the mushy region, say at z = 0.01m, grows when t passes 15s. According to the definition given by eq. (41), this exhibits thermal runaway. The field distribution at t = 24s decays quickly in the water and mushy regions and slowly dies out in the ice region to the right. From the analysis carried out above, one can easily see why it is so difficult to heat any water-containing substance from the frozen state; the latent heat together with the high loss factor of water above the freezing point are in effect a barrier that allows the temperature to rise quickly on one side where T >0°C while suppressing the temperature rise on the other side where T < 0°C. It is noted that the problem of one-dimensional microwave heating of frozen sub stances has also been studied by Coleman [31] using both analytical and numerical approaches. The results he obtained are qualitatively similar to the ones given above. Coleman concentrated in solving only the thermal equation while ignoring the wave equation by assuming exponential field distributions in the water and ice regions. This assumption is acceptable but strictly incorrect because of the field reflection at the interface of the two regions. One interesting result he obtained is that, for some materials with relatively high thermal diffusivity, the mushy region between the thawed and frozen regions can be totally eliminated, especially when the microwave heating is supplemented by convective heating at the surface of the medium. The convergence of the numerical solutions has been tested for all three cases investigated above. For cooked beef juice and Nylon 66 , the numerical solutions converge well even with the length steps and time steps twice as large as the values chosen above but the results are not shown here. Instead, we look at the case of ice-water for which the numerical convergence could be expected to be the worst among the three cases because o f the rapid changes o f many parameters around the R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 58 (x) : h=0.0005m k=0.5s (o): h=0 .001 m k= 1.0s (+) : h=0 .002 m k=2 .0 s (*): h=0.004m k=4.0s T[°C\ h: length step k: time step -10 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 z [m] Fig. 3-9 Convergence o f numerical solution for temperature distribution in water at r = 24s freezing point. Fig. 3-9 shows, nevertheless, a robust convergence of the solution of T for h = 0.001 m and k = 1 s. A further reduction in length and time step yields only a negligible difference. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission . 59 4. TRANSVERSE MICROWAVE HEATING - THERMAL RUNAWAY The one-dimensional microwave heating considered in Chapter 3 is not just o f illustrative value but also it applies approximately to the situations where the dielectric bodies are very large in wavelengths. In such situations, the incident microwave has no other way to reach the interior of the body than penetrating through the surface and dissi pating most of its power there. The power penetration, thermal runaway and all the characteristics of one-dimensional microwave heating depend totally on the dielectric and thermal properties o f the body. In practice, the physical dimensions of most dielec tric bodies subject to microwave heating are of several wavelengths or less. The finite size of the bodies necessitates analysis of the field inside the Fl/ uos in three dimensions. Here, beside the material properties, the body geometry plays an important role because o f resonances or modes that exist in the bodies. The three-dimensional analytical treat ment of the electromagnetic field in a dielectric body is quite complicated even for bodies with regular shapes like circular or rectangular cylir ders. In this chapter, we will again use the wire grid modeling method to compute the power dissipation for several lossy dielectric slabs. Such bodies are chosen to illustrate some frequently observed phenomena like edge heating and thermal runaway of a cheese slice in a microwave oven. Fortunately, _hese phenomena can be explained by surface wave theory without involving the complication of the three-dimensional analysis o f dielectric resonators. In contrast to the case of one-dimensional heating, the thicknesses o f the slabs to be considered a*e just a fraction of one wavelength so that the incident fields can easily penetrate the slabs and non-uniform heating only happens in the transverse planes. We therefore refer to this kind o f heating as transverse microwave heating. 4.1. Surface wave theory The established theory of surface waves on a dielectric slab will be summarised R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 60 here because of its direct applicability to the problem at hand. Consider the infinite flat dielectric slab in free space, shown in Fig. 4-1. Wave propagation is taken to be in the x —v r direction and the x - variation of the fields is assumed to take the form X ( x ) x e , where yx is the complex propagation constant in the x direction. As in the analysis of a parallel-plate waveguide, the wave fields are assumed to be independent of the y coordinate, so that all y - derivatives in M axwell’s equations are zero ( dy = 0 ). With these assumptions, Maxwell’s equations take the following form, in which e = £i inside the slab and £ = £o outside: dHx dEx - f c - + IxH z = -= £ - -jJxHy = ja x E z dHy + = ju \iH y - j y xEy =-ja>\iH2 dEy = j coe£z (4-1) = j($\xHx The permeability p is taken to be the same everywhere, with p = Po- The field inside and outside the slab must have the same propagation constant yx in order that the field con tinuity conditions can be satisfied at all points on the two sides o f the slab suiface. z Eo t I/ E’ y vc.* Propagation ----------- ► \ \ y -t Fig. 4-1 Edge view o f infinite slab supporting surface waves R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 61 Maxwell’s equations (1) are readily separable in two groups, the TE group consisting of the second and the third equation in the right column together with the first equation in the left column, and the TM group consisting of the remaining three equations. 4.1.1. TE waves From M axwell’s equations (1) the relevant TE field components are Hx , Ey and Hz. As in the analysis of rectangular waveguide, we use the separation-of-variables tech nique, assuming that Hx(z,x) = Z(z)X(x). The solution for Hx inside the slab takes the form: Hx l = C {Z (z)e~v (4-2) where C i is an arbitrary constant and Z (z) is proportional to sinh(n \ z ) or cosh(w \ z ) or a linear combination of the two functions. The solution is most conveniently expressed in terms of these funtions because u t , defined by the separation relation u \ = -(co 2p.£1+y^) (4-3) in general, is a complex constant. Even Modes. Choose Z (z) = sinh(u i z). Inside the slab, the complex solution for TE-even modes is Hx i(z,x) = C isinh(w iz)e ^ (4-4) £ y i ( z ,x ) = ^ ^ C 1c o s h ( u 1z )e (4-5) J ti . Hz \(z,x) = — C x C o s h ( u i z)e YxZ “l Outside the slab, for z>t the solution is of the form R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. (4-6) where Co is a constant and uq is defined by uq = -(co : 'i£o+ ^) (4 10) Note that Im («o) > 0 “ a Re {«o ) > 0 to assure a vanishing field for z — At the boundary surface z = r, the condition Hx y = Hxq gives C isin h (u !r) = Co^ u ' (4-11) The boundary condition Ey \ = Ey o at z = r gives Ci ui cosh(uiz) = Co - Ur,i e uo (4-12) Eliminating C \ I C q gives witanh(uir)+uo = 0 (4-13) With u\ and u q given in eqs. (3) and (10) respectively, the transcendental equation (13) can be solved to give the propagation constant yx for a given mode. The Muller method [25] is most convenient for seeking a complex root of the equation. Odd Modes. Choose Z(z) = cosh(u i). The soludon for TE-odd modes inside the slab is Hx \ (z,x) = Ci Cosh(ui z)e' yiX Ey i (z,x)= ] sinh(«\z )e ^ llX u! R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission. (4-14) (4-15) Outside the slab, the solution is given by eqs. (7)-(9) as before. Matching the tangential fields at z = t gives the transcendental equation for TE-odd modes u 1coth(t' it>fUo = 0 (4-17) It turns out that, for thin slabs with thickness much smaller than the intrinsic wavelength of the dielectric medium, only the lowest even mode TEo can propagate. 4.1.2. 7 M M/.<7v*».r From M axwell’s equations (1), the relevant TM field components are Ex , Hy and Ez. Following the derivation procedure as in the TE case, we get the solutions with asso ciated transcendental equations for TM even and odd modes. Even M -ie s. The complex solution forTM -even modes inside the slab is Ex \(z,x) = C \S\n\i{u\z)e v Hy i ( z , x ) = J (0£l _yr C;Cosh(u!z)e f (4-18) (4-19) “1 £, i( z ,x ) = —— CiCosh(uiz)e YrJ: Ul (4-20) where yx is the root or the transcendental equation Ui £i and u \ , uq tanh(u 11)H Uq =0 Eo are defined as above. Odd Modes. The complex solution for TM odd modes inside the slab is R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. (4-21) 64 £ ^ 1(z,x) = CiCOsh.''w12)e Yp: (4-22) (4-23) Ez i(z,x) = — ^ -C 1sinh(u 1z )f ^ Ml (4-24) where yx is the root of the transcendental equation ill (4-25) El and u i , Uq are defined as above. 4.2. Edge heating in infinitely long dielectric strips For lossless dielectric slabs, yx is imaginary and the surface wave propagates without attenuation. For lossy slabs, yx is generally a complex constant and the surface wave dies out with the attenuation constant a x, where cc* = Re {yx }. Consider the infinitely long dielectric strip shown in Fig. 4-2. The strip can be thought o f as a segment of the infinite slab shown in Fig. 4-1. Assume the field in the strip can be expressed in terms of the surface wave modes in an infinite slab. We are interested in investigating the field variation inside the strip when a wave is incident upon it. It is o f interest to con sider the two cases, TE and TM: 4.2.1. The TE case Let the strip be illuminated by some incident wave with the electric vector field parallel to the edges of the strip. The tuiai field inside the strip consists o f the "forced wavv" and surface waves. By "forced wave", we mean the wave that has a variation only in the z direction b' i no variation along the x-direction. The existence o f this wave can be verified by modifying e -YrX with e ^ + C, where C is a constant, in eqs. (4)-(6) and R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 65 y -a Fig. 4-2 Geometry o f infinite dielectric strip with TE surface waves Fig. 4-3 Normalized power density distribution in infinite cheese strip (er = 5 4 -/2 1 ) due to TEq surface waves. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 66 still satisfying the Maxwell’s equation given in eq. (1). The other waves are TE surface waves which are induced by diffraction of the incident wave at the edges of the strip. (Because of the infinite length of the strip, the total E-field has no x — component and hence no TM wave is present in the total fieldV The TE surface waves propagate either in the positive or in the negative x direction and the propagation constants for these waves are determined by eq. (13). Usually, only the lowest even modes TEo can propagate on a thin dielectric strip; other non-propagating modes may contribute to the field at the edges but they most probably die out very quickly. Richmond [26] investigated the same prob lem using Galerkin’s method to determine the amplitudes of different modes existing in the strip. He found that it is accurate to express the field in thin strips with TEo waves and the forced wave only. If the forced wave is small compared to the TEo waves, the field variation along the x direction can be approximately expressed by X (x) = A e ^ + B e (4-26) where A and B are constants. For symmetric excitation, A = B and we have X ( x ) = 2A c o s h ^ x ) (4-27) The normalized power density distribution is given by IX CO | 2 2 |cosh(yxs ) | t = '— ; |T |X (n ) | 2 |cosh(Yra ) | (4-28) Obviously, if ct* is large, that is a xa :» 1, the field attenuates quickly from the edges so the field magnitudes at the edges will be much larger than at the center of the strip, i.e. w*. have edge heating. Otherwise, if a xa < 1, p(x) may have a large local maximum at the center, i.e. we have center heating. Example. Dielectric properties of cheese at / = 2.45GHz and T —50°C are e / = 54 and ?/■’ = 21. For a cheese slab of thickness 2r = 2mm, the propagation constant for the TEo surface wave is found from eq. (13) to be yx = 44 + j 139 m-1. A plot of p (x) for a cheese R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 67 strip of width 10cm is ^ ° w n in Fig. 4-3. It is seen that the power density dies out quickly from the edges; at distances 1cm from the edges, the power density is reduced to less than half of the maximum value. It is worth noting that the nature of decaying surface waves on a thin lossy strip agrees well with the edge behaviour of surface current on a strip conductor. It can be shown that for a perfectly conducting strip, the attenuation constant for TE surface waves is infinitely large, which means the surface current dies out "instantly" in the vicinity of the edges. Since the surface current at the center of the strip is nonzero, the surface currents at the edges are forced to infinity. The singularity of edge current is well known and was observed by several authors [27], [28] when they used different methods to cal culate the surface current distribution on a perfectly conducting strip being illuminated by plane waves with their electric field vectors parallel to the edges of the strip. \ 4.2.2. The T M case ; Consider the same dielectric strip as before; the electric field vector is now perpen dicular to the edges of the strip. It is clear that only TM surface waves are induced in this ! case. For strips with high dielectric constant or loss factor, it is found that the propaga tion constant fc. TM waves is very close to that o f plane wave in free space. For the ; cheese strip in the above example, the propagation constant for TMo wave is found from eq. (17) to be yx = 0.0038 + j 51.42 m_1. For comparision, the propagation constant of a plane wave in free space ;s y’Po = yco^poEo = j 51.35 m_1 at / = 2.45GHz. In contrast to the TE case, the attenuation constant for the 1 ?vl wave is very small, which does not lead to edge heating as in TE case. For normal illumination, the two TM waves propagating in the positive and negative x directions have the same magnitude arid the resultant field distribution has a form similar to the standing wave pattern on an electric dipole in free r space. r R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. ] i 68 4.3. Edge heating in rectangular dielectric slabs The infinite strip considered above was useful to illustrate the edge heating phenomenon in lossy dielectric. In practice, any dielectric material subject to microwave heating has finite dimensions and the heating pattern for such a three-dimensional object becomes understandably more complicated. It is most natural in the study of edge heating to consider a rectangular dielectric slab as shown in Fig. 4-4. The slab can be thought of as being a segment of the strip considered above. However, no straightforward derivation can be made for the slab from the analysis of the strip because the assumption dy = 0 is no longer valid. An appropriate analytical treatment of the field inside the slab would involve the theory o f dielectric resonators which is beyond the scope of this thesis. Even if the propagation constants for different modes in the slab cavity could be determined analytically, the simultaneous existence of many modes with unknown magnitudes would not help one in gaining an insight into the true field existing in the slab. In this section, we will instead use the wire grid modeling method to study the heating pattern for some dielectric slabs of food materials with known dielectric properties. The surface Z Fig. 4-4 Geometry of rectangular dielectric slab R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 1 69 wave theory can still be used approximately to choose the approptiate dimensions for the slabs and to discuss the results obtained. To demonstrate edge heating in a rectangular slab, we choose to model three slabs with the same dimensions but different food materials. The dielectric properties for these materials are estimated from the data given by Ohlsson and Bengtsson [29]. The related propagation constants and wavelengths for these materials are listed in Table 4-1, where _ y = ct + v'p = (-cd2 [Io£o)1/2 is the propagation constant of a plane wave in an infinite dielectric medium _ Is = a j + is the propagation constant of a TEo surface wave in an infinite dielectric slab of thickness 2 1 = 2 mm It should be noted that the listed values of e / and er" are just estimates at / = 2.45GHz and T = 50°C for the chosen materials because of the limited quantity of available data. This, however, does not affect the theoretical aspect of the study. The phase-shift con stants P are of the same order for the three kinds of food but there are significant differ ences between the attenuation constants for these niaterials, with ham being the lossiest one. Similar observations can be made for P5 and a*. The choice of the materials has in fact been based on these criteria to achieve slabs with the same physical dimensions as well as the same wire grid model (but different loads). The results obtained for these slabs would then be most easily compared and interpreted. Since the surface wavelengths Xs of all three materials are about 0.045 m, an appropriate length for the slab would be 2a = 0.064 m which is about 3/2 Xs. This is certainly the minimum length that would allow us to see the edge heating phenomenon in the slab. The width 2b of the slab can be chosen more freely but it should be much larger than the thickness o f the slab R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 70 to keep the assumption of surface waves propagating on thin slab valid. We choose b = a/2 = 0.032 m. The thickness of the slab has been chosen to be 2 mm, which is about 1/8X or the maximum allowpble wire grid spacing. The dimensions of the full wire grid model for the slab are then (2ax2£x2c) = (32x16x1) d, where d = 2 mm is the grid spac ing. Using symmetry, the size of the quarter wire grid is (16x8x1) d. The source dipole is y-oriented and placed lm above the slab. This large distance together with the utilized ME-symmetry guarantee a virtually uniform plane wave incident normally upon the slab. The whole configuration is shown in Fig. 4-5. The wire grid model for each material is loaded with corresponding lumped impedances calculated from eqs. (2-4)- (2-6). Volume resistance Rv and capacitance Cv are listed in Table 4-1. The computed power density distributions for the three cases are shown in Figs. 4-6. The expected edge heating is clearly shown in Fig. 4-6a, where the power density in the cheese slab has its maxima at the two edges (x = ±a) and decays toward a minimum at the center (x = 0). The yvariation of the power density takes the form of a standing wave with its maximum at y = 0. It is seen here that the concept of decaying surface waves is still applicable to ex p iJ^ ; edge heating phenomenon although the attenuation and phase-shift constants in the x direction are different from those in case of the infinite strip. The evidence for surface waves is more obvious in the less lossy pineapple slab for which the power den sity is shown in Fig. 4-6b. Here, we again see two maxima at the edges but what distin guishes this case from the previous one is the local maximum at the center of the slab. It should be noted from the previous section that the field distribution on an infinite strip always has a local maximum at the center if the strip is illuminated symmetrically. How ever, if the strip is very lossy or wide, the surface waves will effectively die out at the center and the local maximum becomes negligible. The same can reasonably be stated for a rectangular slab; we see that no local maximum is present at the center of the cheese slab because the surface wave’s attenuation constant a* for cheese is much larger than that for pineapple (see Table 4-1). The distribudon of power density in a rectangular slab R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. 71 Table 4-1 Dielectric properties for some food materials at / = 2.45 GHz and related propagation constants, wavelengths and volume lumped loads for wire g ud model Food £r Y [m-1 ] X [m] Yj [m l ] Xj [m] Rv [« ] Cv [pF] Cheese 54-j21 72+J384 0.0164 44+jl39 0.0450 88 1.88 Pineapple 51-j9 32+j368 0.0171 19+jl45 0.0434 204 1.80 Ham 50-j38 130+j386 0.0163 81+jl31 0.0470 48 1.76 (0 .02, 0 . 02, 1.00 ) 2 cm source dipole Fig. 4-5 Quarter grid model for (64x32x2) mm dielectric slab with source dipole R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. 72 a) b) Fig. 4-6 Power density distribution in a) cheese slab b) pineapple slab c) ham slab o f (64x32x2)mm R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 73 ’s admittedly much more complicates than in the case of a corresponding infinite strip because it is a combination o f many modes that exist in the slab. With ham being the lossiest dielectric considered, the power density distribution in Fig. 4-6c shows an extreme case of edge heating where the power density highly concentrates at two 15mm wide strips at the edges. As a final investigation on edge heating, consider a square cheese slab of (16x16x2) d being illuminated by a uniform plane wave. The whole configuration shown in Fig. 4-7 is identical to the one shown in Fig. 4-5 except the length of the slab has been reduced to 2a = 16 d = 32 mm which is only 3/4 Xs . The power density distribu tion is shown in Fig. 4-8. In contrast to the cases studied above, the power distribution has its maximum at the center of the slab, which can be reasonably explained in the same way as before. Due to the relatively small length of he slab, the surface waves travelling from the edges have not decayed significantly before they interfere and build up the standing wave peak at the center. One special characteristic of power density in the slab is the peaks that are highly localized at the comers. Simulations of microwave heating in other slabs showed that the peaks only show up for small slabs but not for those whose dimensions are in order of wavelengths (the results are not shown). Similar observations have been made by Ko and Mittra [28] for surface currents *n perfectly conducting plates. There seems to be no exact explanation for this phenomenon at the present time. 4.4. Thermal runaway in transverse microwave heating Just as in the case of one-dimensional heating, initial non-uniform heating can lead to thermal runaway in the material being illuminated transversely. In microwave transverse heating, the source of initial non-uniform heating can be either uneven illumi nation or an edge heating effect. Whatever the reason, the temperature dependence o f the dielectric properties determines the tendency of a material to develop thermal runaway. Thermal runaway is typically observed by heating a rectangular cheese slice in a microwave oven. Usually, there are two parallel edges being overheated. As time passes, R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 74 (0.02 , 0 .02, 1.00) 2 cm source dipole d = 2 mm y Fig. 4-7 Quarter grid modei for (32x32x2) mm cheese slab with source dipole f i a Fig. 4-8 Power density distribution in (32x32x2) m m cheese slab R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. 75 the borders between these edges and the rest become more and more well defined; the edges finally melt and degrade while the rest is bareiy warmed up. It is o f interest to con sider the cheese slab with the wire grid model shown in Fig. 4-5. The temperature depen dence of the dielectric properties of cheese can be described by the linear approximations [291: e / = 54 + ^ ( 7 - 5 0 ) er" = 21 + - j y ( r - 5 0 ) in the temperature range T = 50 —» 125°C. After an initial uneven heating, the slab becomes dielectrically inhomogeneous due to the non-uniform temperature distribution in the slab. The wire grid model for the slab must therefore be loaded non-uniformly. The new lumped impedances for the wire grid model can be calculated as described previously in Chapter 2 with the assumption that the temperature and hence the dielectric properties are constant over any given elemen tary dielectric cube of the model. With the new loads, a new power distribution must be computed. With the assumption that the heat exchange in the slab is negligible due to the quick heating process, the dynamic heating pattern in the slab can be simulated as fol lows: 1/ Set up the wire grid model for the homogeneous slab. 2/ Compute the power density distribution 3/ Compute the temperature rise in the slab AT = pC Usually, pC is constant with temperature and the time step can be scaled to give a maximum temperature rise ( A = — • The temperature rise is then generally Pc given A T = - 2 - ( A T ) max. P m ax R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. by 75 the borders between these edges and the rest become more arid more well defined; the edges finally melt and degrade while the rest is barely wanned up. It is of interest to con sider the cheese slab with the wire grid model shown in Fig. 4-5. The temperature depen dence of the dielectric propenies of cheese can be described by the linear approximations [291: e / = 54 4 - ^ ( 1 - 5 0 ) e / ' = 21 + - ^ - ( r - 5 0 ) in the temperature range T = 50 —» 125°C. After an initial uneven heating, the slab becomes dielectrically inhomogeneous due to the non-uniform temperature distribution in the slab. The wire grid model for the slab must therefore be loaded non-uniformly. The new lumped impedances for the wire grid model can be calculated as described previously in Chapter 2 with the assumption that the temperature ana hence the dielectric properties are constat over any given elemen tary dielectric cube of the model. With the new loads, a new power distribution must be computed. With the assumption that the heat exchange in the slab is negligible due ro the quick heating process, the dynamic heating pattern in the slab can be simulated as fol lows: 1/ Set up the wire grid model for the homogeneous slab. 2/ Compute the power density distribution p A' 3/ Compute the temperature rise in the dab .jT --—— ' Is iallv, ->C is constant DC with temperature and the time step can b" 'T e d to g ve a m'ximum t mperature CmaxAr rise (AT)mM = ----- -— . The temper:-1 ' pC rise is d u n . lerall • given by A7 = - £ - ( A r w P nax R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission. 76 4 / Calculate the new distribution of dielectric properties and hence the new lumped loads for the model. 5/ Go to 2/ for the next time step. If the time step is small such that the temperature rise and hence the dielectric change are small then the simuladon should approach the rea* evolution of power density in the slab. For the cheese slab, a maximum temperature rise is chosen to be 25°C for each time step. Assuming yC = 1.0-10'J/(m3oC) the results are shown in Fig. 4-9. Thermal runaway is evident by the rising power density at the edges and the remaining low power density it the center. Thermal runaway is expected in :his case because of the rising loss factor er" with temperature, which means the hot regions in the slab get lossier and become hotter and so on. Furthermore, the power carried by the <urface waves travelling from the edges to the center gets absorbed more and more n r r the edges leaving virtually nothing for the center. However, thermal runaway is not very strong in this case probably because the dielectric constant e / is also rising with tempera ture, which in effect reduces the penetrability of waves into the hottest regions. Admit tedly, the investigation conducted above is very' sirr~'e and is only one of many cases that show thermal runaway; a rigorous treatment of thermal runaway would require much more research in this virtually unexplored area. R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. 77 a) a b) c) Fig. 4-9 Evolution of power density distribution in cheese slab: a) r = Os, b) r = 76s, c) t — 146s R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. 78 CONCLUSIONS The wire grid modeling method has been shown to be successful in modeling con ductors and dielectric materials. Simulation of microwave heating of a dielectric in a small, deep rectangular waveguide model shows that the numerical results obtained by the method are comparable to know: analytical solutions. The method was then used to investigate the microwave heating properties of various dielectric models, including edge heating and thermal runaway phenomena in common food materials. The predicted heat ing patterns in these models are consistent with the practical heating patterns of food conanonly observed in commercial microwave ovens. General aspects of one dimen sional microwave heating have also been studied to give some insights into the physics of the thermal runaway effect. In particular, the finite difference method was used to solve coupled wave-thermal equations witti some simulations of one-dimensional microwave heating revealing the tendency of several dielectric materials to develop thermal runawav With the ability of modeling one- two- and three-uimensionai objects for a wide variety of materials, the wire grid modeling method is shown to be a powerful tool, not only for the research in microwave heating but also for solving general electromagnetic scattering problems. Despite the limitation that only small objects of several wavelengths can be modeled with the wire grid modeling method, it is hoped that the rapid evolution o f the modem computer technology will enaLe larger models to be run on machines with higher speed and larger memory and this will certainly open up a new spectrum of prob lems that can be solved by the method. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 79 APPENDIX A IM PLEM ENTATION O F ELEC TR IC AND M AGNETIC SYM METRIES IN TH E TH IN -W IRE PROGRAM In electromagnetic analysis, the theoretical perfect electric conductor is frequently used to get the simple boundary condition nxE = 0 or E tan = 0. As a dual to the perfect electric conductor, the hypothetical perfect magnetic conductor is the one that satisfies the boundary condition -n x H = 0 or H tan = 0. The images of electric current elements normal and parallel to these perfectly conducting planes are shown in Fig. A -l The concept of electric and magnetic symmetries is extremely useful in the perfect electric conducting plane (E-planc) positive image negative image pcncct magnetic conducting plane (M-planc) negative image positive image Fig. A -l R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited without p erm ission . 80 reduction of unknown dipole current modes when solving a problem using the thin wire program. If the wire structure is symmetric with respect to the X=0 and Y -0 planes for instance, electric and magnetic symmetries can be imposed on these planes to reduce the unknown currents to the ones flowing in the first quadrant only. The currents in other quadrants are images of the first quadrant currents. The physical interpretation o f impos ing a symmetry on a plane is that the wire structure is symmetrically excited or illuminated by electromagnetic sources around that plane. Consider the simple wire structure consisting of two dipoles shown in Fig. A-2a. System equations for the wire dipole 1 i~ dipole 2 i—. i i V'i V2 a i i Fig. A-2a i l — >dipole E-planc -----------------------------------------image Fig. A-2b -T, -c------- A dipole V'. M-planc — — — — — — — — image Fig. A-2c c :::::: -----I R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. 81 structure are Z i i / i + Z i 2/ 2 = Vi (A-l) ^ 21^2 "h ^ 22^2 = ^ 2 (A-2) We always have Z 12 = Z 2 i. And Z n = Z 22 if the dipoles are identical. V l and V 2 are excitation voltages. Electric symmetry consideration. If V x = - V 2 then it readily follows 1 1 = - / 2 and the configuration corresponds to the one shown in Fig. A-2b with If the dipole and its image coincide, that is if the dipole is on the symmetry plane then Z\ \ = Z 12 and I \ is undetermined. In general, therefore, when solving for a thin wire problem using electric symmetries, the segments lying on these planes must be removed to avoid the above problem. This is possible because the total currents on these segments are zeros and removing them does not effect the solution for the other currents on the wire structrure. Another advantage is that the number of unknown currents will be then reduced. Magnetic symmetry consideration. U V X = V 2 then we have I x =1 2, which corresponds to the configuration with mag netic symmetry shown in Fig. A-2c, and Vi '1 = 7 — 7 Z 11 + -M2 Now, if the dipole is on the symmetry plane then Z \ \ = Z i 2 and R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. (A-4) 82 (A-5) 2Z 11 2 where I \ is the current flowing in the dipole in free space. If the dipole is loaded with Z i , / ^1— 2Zn + Z l 1 (A- 6 ) and the dipole in free space is loaded with Z' \ , Vi /' i = (A-7) Z n +Z\ Since the total current on the segment lying on the symmetry plane is 21 [ which should be equal to l \ for all values of V i , it follows that Z t = 2Z 'i. When using symmetries, this means the segments lying on magnetic symmetry planes should be loaded witn twice the impedance these segments are loaded with when no symmetry is used. For dipoles perpendicular to symmetry planes, it is worth noting that currents are able to flow normally through an electric plane but are shorted in the case of a magnetic plane. In practice, this means, if a dipole is placed perpendicular next to ar. electric plane, an extra current mode must be set up at the plane. In contrast, no extra current mode is needed if the plane is magnetic. See Fig. A-3. M-plane E-plane r> 11' II ' u n extra current mode1 n II “U r / " , ly Fig. A-3 R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 83 APPENDIX B COM BINATION OF SYM M ETRIES The use of electric and magnetic symmetries for a symmetrical wire structure implies that it must be symmetrically excited. However, nonsymmetrical excitation can be achieved by combining different symmetry configurations together since electromag netic quantities are linear and additive. Consider the simple wire structure shown in Fig. B-la. Dipole 1 is fed with current I and dipole 2 with a zero current. The configuration is equivalent to the linear combination of the electric and magnetic symmetries as shown in Fig. B -lb. Combination of two-dimensional symmetries can also be carried out and is shown in the following. The four basic two-dimensional symmetry configurations with x- and y-oriented current elements are shown in Fig. B-2. Let x ,, x 2, x 2 and X 4 be the coefficients of the basic configurations MEx, EMx, MMx and EEx respectively in the combination to create a desired excitation configuration with x-oriented currents having the values lclt k 2, £3 and k.4 in the 4 quadrants. Referring to Fig. B-2, we can write the combination equations in matrix notation as follows, / , =1 Fig. B -la / 2 =0 1/2 1/2 E-plane + a M-plane Fig. B-lb c R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 84 M M M Ex MEy EM x EMy M M M MMx MMy — M M EEx EEy Fig. B-2 R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 85 1 1 -1 il 1 1 1 1 -1 -11 -1 - 1 1 - 11 ' 'k{ *1 ki *2 “* (B -l) *3 *3 X4 k4 Solving for x x to x 4, we get ’ 1/4 -1/4 1/4 —1/4 *2 1/4 1/4 1/4 1/4 “ 1/4 -1/4 -1 /4 1/4 *3 1/4 1/4 -1 /4 -1 /4 x4 ' k f ki ki (B-2) k4 Similarly, let y i , >’2 , >3 and y 4 be the coefficients for MEy, EMy, MMy and EEy configurations respectively. And let l \ , I 2 , Ij, and l4 be the values of y-oriented currents in the 4 quadrants. Referring to Fig. B-2, we have the matrix equation 1 1 1 -1 -1 il 'yi 1 yi 12 1 y h 1 ya 3 • / r (B-3) U which gives yi yi y3 y4 ' 1/4 1/4 1/4 1/4 1/4 -1/4 1/4 -1/4 1/4 1/4 -1/4 -1/4 1/4-1/4-1/4 1/4 -li h h U R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. ('3-4, 86 REFERENCES [1] W.L. Stutzman and G.A. Thiele, "Antenna theory and design", John Wiley & Sons, New York, 1981, pp. 356-364. [2] J.H. Richmond, "Radiation and scattering by thin-wire structures in the complex frequency domain", National Technical Information Service, Springfield, VA 22151, NASA CR-2396, May 1974. [3] C. Aversa, "Microwave properties of wire grids", B.A.Sc., thesis, Div. of Engineering Science, University of Toronto, Toronto, Canada, 1988. [4] A.C. Metaxas and R.J. Meredith, "Industrial microwave heating", Peter Peregrir.us, London, 1983, pp. 5-24. [5] N.E. Bengston and P.D. Risman, "Dielectric properties of foods at 2.8 GHz as determined by a cavity pertubation technique", J. Microwave Power, vol. 6 , no. 2, p. 107, 1971. [6 ] E.C. To, R.E. Mudgett, D.I.C. Wang, S.A. Goldblith and R.V. Decareau, "Dielec tric properties of cood materials", J. Microwave Power, vol. 9, no. 4, p. 303, 1974. [7] R.F. Harrington, "Field computation by moment methcds", Macmillan, NewYork, 1968. [8] R. Mittra, ed., "Computer techniques for electromagnetics", Pergamon Press, Oxford, 1973. [9] K.S.H. Lee, L. Marin and J.P. Castillo, "Limitation of wire grid modeling o f a closed surface", IEEE Trans. Electromag. Compat., vol. EMC-18, no. 3, pp. 123129, Aug. 1976. [10] J. Moore and R. Pizer, "Moment methods in electromagnetics", Research Studies Tress, New York, 1984, secs. 1.1.2 and 6 4. [11] A.C. Ludwig, "Wire grid modeling of surfaces", IEEE Trans. Antennas Propagat., vol. AP-35, no. 9, pp. 1045-i048, Sept. 1987. [12] E.C. Jordan and K.G. Balmain, "Electromagnetic waves and radiating systems", * entice-Hall, Toronto, 1968, pp. 470-473. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission. 87 [13] M.A. Tilston, "Thin-wire reciprocal multiradius implementation o f the elec tromagnetic moment method", Ph.D, thesis, University o f Toronto, Toronto, Canada, 1989. [14] P.C. Clemmow, "The plane wave spectrum representation o f the electromagnetic field", Pergamon Press, Oxford, 1966. [15] D.T. Paris and F.K. Hurd, "Basic electromagnetic theory", McGraw-Hill, 1969, New York, pp. 367-370. [16] J.R. Wait, "Electromagnetic waves in stratified media", Pergamon Press, New York, 1962, pp. 88-91. [17] in [16], pp. 100-105. [18] D.M. Young and R.T. Gregory, "A survey of numerical mathematics", AddisonWesley Publishing Company, New York, 1972, pp. 473-477. [19] "Numerical recipes: the art of scientific computing", Cambridge University Press, 1986. [20] in [18], pp. 357-358. [2 1 ] in [18], pp. 1086-1088. [22] in [18], pp. 586-587. [23] in[4], p. 55. [24] E.N. Ganic, J.P. Hartnett and W.M. Rohsenow, "Handbook of heat transfer funda mentals", McGraw-Hill, New York, 1985. [25] D.E. Muller, "A method for solving algebraic equations using an automatic com puter", Mathematical Tables and Aids to Computation, vol. 10, 1956, pp. 208-215. [26] J.H. Richmond, "Scattering by thin dielectric strips", IEEE Trans. Antennas Propagat., vol. AP-'J>, no. 1, op. 64-68, Jan. 1985. [27] S.M. Rao, D.R. V/ilton and A.W. Glisson, "Electromagnetic scattering by surfaces o f arbitraiy shape , IEEE Trans. Antennas Propagat., vol. AP-30, no. 3, pp 416- R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission. 88 417, May 1982. [28] W.L. Ko and R. Mittra, "A new approach based on a combination of integral equation and asymptotic techniques for solving electromagnetic scattering prob lems", EEEE Trans. Antennas PropagaL, vol. AP-25, no. 2, pp. 187-197, March 1977. [29] Th. Ohlsson and N.E. Bengtsson, "Dielectric food data for microwave sterilization processing", J. Microwave Power, vol. 10, no. 1, 1975. [30] J.T. Mayhan, "Characteristic modes and wire grid modeling", IEEE Trans. Anten nas PropagaL, vol. 38, no. 4, p. 459, April 199G. [31] C.J. Coleman, "The microwave heating of frozen substances", Appl. Math. Model ling, vol. 14, pp. 439-443, Aug. 1990. R ep ro d u ced with p erm ission o f the copyright ow ner. Further reproduction prohibited w ithout p erm ission.

1/--страниц