Numerical Modeling of Microwave Heating

The present study compares the temperature distribution within cylindrical samples heated in microwave furnace with those achieved in radiatively-heated (conventional) furnace. Using a two-dimensional finite difference approach the thermal profiles were simulated for cylinders of varying radii (0.65, 6.5, and 65 cm) and physical properties. The influence of susceptor-assisted microwave heating was also modeled for the same. The simulation results reveal differences in the heating behavior of samples in microwaves. The efficacy of microwave heating depends on the sample size and its thermal conductivity.


Introduction
Microwaves are now well-recognized mode of heating material.In the beginning (late 40's) 2.45 GHz microwaves were first used for heating food products up to 100-150°C [1].Subsequently, they were used in rubber industry for aiding vulcanization [1].Later on, microwaves were utilized for heating ceramic bodies [2].It was reported that microwaves interact both with bulk as well as particulate ceramics and can cause rapid increase in their temperature to up to 1700°C [3].The coupling of particulate ceramics elicited keen interest among researchers as it provided means to rapidly sinter the ceramic powder compacts [4] without significant microstructural coarsening.Unlike conventional heating, microwaves directly couple with the material, which leads to rapid and uniform heating.In the context of sintering this is particularly important.Over the years, a variety of materials such as carbides, nitrides, complex oxides, silicides, and apatites, etc. have been synthesized using microwaves.Despite the widespread application of microwaves for thermal processing of materials, their use till recently restricted for metals.This was on account of the fact that bulk metals tend to reflect most of the microwaves incident on them, and only heat over a localized surface.Such an interaction frequently results in arcing that eventually causes the failure of the magnetron.However, in 1999, for the first time Roy et al. [5] demonstrated that most of the particulate metals very effectively couple with microwaves and can be heated to high temperatures at significantly faster rates.This has opened the possibility of consolidating particulate metal compacts to full density through microwave sintering [6].Several researchers [7][8][9][10][11][12][13][14][15][16] in recent years have demonstrated the efficacy of microwave sintering in a range of materials that include Cu, Sn, bronzes, stainless steel and tungsten heavy alloys.
Despite the widespread use of microwaves for heating materials, there is still no unanimity regarding its interaction mechanism with the materials, particularly, metallic ones [17][18][19][20][21].The main loss mechanisms during microwave-material interaction are electric, conduction (eddy current), hysteresis and resonance (domain wall and electron spin (FMR)).It is often difficult to ascertain which loss mechanism, or combination of mechanisms is occurring for a particular sample in given conditions.The different mechanisms do however have different dependencies on certain properties such as sample type and microstructure, frequency and temperature.
Recently, Mishra et al. [22] modeled the heating of particulate metal compacts by solving the Maxwell equation and using a 2D-Finite Difference Time Domain (FDTD) model.While their model showed good correlation between the experimental and simulated temperature profiles of Cu, Sn and W, it had several limitations.Their model ignored the effect of compact size and its thermal conductivity and assumed a uniform temperature distribution within the compact.Furthermore, it failed to compare the thermal profile of microwave heated compact with those achieved using a conventional (radiatively heated) furnace.In case of microwave heating, since the sample per se acts as the source of heat, hence, it is hypothesized that the compact should have an inverse thermal profile with the central (inside) portion having higher temperature than the peripheral regions.Furthermore, at times, a low temperature the coupling between microwaves and the compact is poor.To address these issues, in practice, microwave heating is usually done by surrounding the compact with susceptors such as SiC that gets preferentially heated at lower temperatures.
The present study, therefore systematically investigates the evolution of temperature profile with time within a cylindrical material as a function of its size (0.65 cm, 6.5 cm, 65 cm) and conductivity under both conventional as well as microwave heating mode.The latter has been examined both in the presence and absence of susceptor.The objective of current study is to device a well defined strategy for the selection of most effective sintering process (conventional versus microwave) for the given conditions.To achieve this, 2-D finite difference mathematical model has been formulated to examine the role of various process variables.

Mathematical Formulation
For the present simulation, the material was considered to interact with 2.45 GHz multimode microwaves.Fig. 1 shows the schematic of a sample configuration in a multimode microwave furnace.For facilitating the modeling an axio-symmetric (cylindrical) sample geometry was assumed with the radius (R) to height (H) ratio of unity.Three different radii (0.65, 6.5 and 65 cm) cylinders were considered in the present study for simulation.To simulate the effect of thermal conductivity on the heating behavior for each size the simulation was performed for a low and high thermal conductivity condition, respectively.
The present configuration was constructed so as to enable the placement of SiC susceptors for assisting microwave heating.Microwaves penetrate and propagate through a dielectric material like SiC.This generates an internal electric field (E) within a specific volume, which in turn, induces polarization and movement of charges.The resistance to these induced motions due to internal, elastic and frictional forces attenuates the electric field.These losses result in volumetric heating.The resulting electromagnetic power absorbed per unit volume (P EM ) [Wm -3 ] by a material is given by [23]: where, ω is frequency of the microwave radiation (2.45 GHz) ε 0 is absolute permittivity of free space (8.85×10 -12 C 2 N -1 m -2 ) ε r ′′ is relative imaginary component of the dielectric constant µ 0 is absolute permeability of free space (4π×10 -7 NA -2 ) µ r ′′ is relative imaginary component of permeability constant E rms and H rms are root mean square values of electric and magnetic field amplitudes, respectively.This heat generation rate per unit volume in the material during microwave sintering can be designated as the source term, S. For the sake of simplification, for constant power and frequency input, the magnetic permeability, electrical conductivity and permittivity are assumed to be linear function of temperature.Elsewhere too, this approximation has been considered [24].The entire source term 'S' therefore can be expressed as: S = A SOURCE + B SOURCE ×T (2) where, A SOURCE and B SOURCE are the constants and are specific to the materials and microwave operating parameters (frequency and power).
The differences between conventional and microwave heating can be described mathematically using unsteady state heat transfer equation with a volumetrically distributed heat source and is expressed as: where, S is local density of heat sources.For conventional heating 'S' is non-existent, whereas for microwave processing it can be determined using the Equations 1 and 2. Having calculated the source term, the heat transfer for the cylindrical geometry was done using a 2D finite difference explicit model.Fig. 2 shows the flowchart of the entire modeling scheme.For formulating the model, the following assumptions have been made (i) A cylindrical coordinate system has been assumed.(ii) Physical properties, such as density, specific heat, thermal conductivity, resistivity are assumed to be uniform throughout the calculation domain.

(iii)
Heat source is uniformly distributed.Rearranging the governing heat transfer Equation 3: For the finite difference modeling, the cylindrical sample geometry is subdivided into radial and axial grids.The total number of grid points in the radial and axial directions has been designated as M and N, respectively (Fig. 3).The grid spacing in both the directions, Δr and Δz, can be expressed as following:  4) can be expressed as: After rearranging, the expression becomes The above equation is valid only for the inner nodes {m: 1, 2, 3 …..M-1 and n: 2…..N-1}.For writing the finite difference equation for the nodes lying on the outer surface and edges the equations have to be reformulated.

Heat Transfer at Side Surface Boundary
In Fig. 4a, the heat balance in control volume centered at location corresponding to node 'O' can be expressed as the summation of the rate of heat flow from various adjacent nodes as following: Rate of heat accumulation in control volume at the side surface ⇒ Rate of heat flow through conduction from {node 1→ node O + 2→O + 4→O + 5→O + 6→O} + Rate of heat flow through convection from 3→O + Rate of heat flow through radiation from 3→O + Rate of heat flow through susceptor 3→O + Rate of heat generation in control volume (11) The above expression can be expressed mathematically as following: Rearranging Equation 12, ( ) ( )

Heat Transfer at Top/Bottom Edge Boundary
As done with the surface nodes, the heat balance in control volume centered at the top and the bottom edges (Fig. 4b) location corresponding to node 'O' can be expressed as the summation of the rate of heat flow from various adjacent nodes to it as following: Rate of heat accumulation in control volume at top/bottom edge ⇒ Rate of heat flow through conduction from {1→O + 2→O + 4→O + 6→O} + Rate of heat flow through convection from {3→O + 5→O} + Rate of heat flow through radiation from {3→O + 5→O} + Rate of heat flow through susceptor {3→O + 5→O} + Rate of heat generation in control volume (14) The above expression is mathematically expressed as: ( ) The above expression can be mathematically formulated as: Rearranging the above, ( )

Error Analysis and Establishing Stability Criteria
The finite difference model involves numerically solving the equation using Taylor series expansion and hence consists of some inherent errors [25].One such error is the truncation error which can be attributed to the computation of derivatives that are used in place of the finite difference equations.Another prevalent error in such formulation is the numerical round-off error that is a consequence of finite significant digit restriction [26].For the explicit forward Euler method used in the present study, it is therefore pertinent to establish the stability criterion.The time step, Δt, for the temporal temperature simulation within the sample is expressed as [27] For the present simulation, assuming a 50×50 (M,N) mesh, irrespective of the sample size and its conductivity Δt ≤ 0.0015s was found to satisfy the stability criterion.

Calculation Procedure
Tab.I provides a summary of the physical properties and other parameters used as inputs for the simulation.To obtain the thermal profiles within the compact at various time intervals, the mathematical formulation of the heat transfer model was done using a FORTRAN based computer program run under UNIX environment.The simulation was run for 180 seconds time duration for both conventional as well as microwave heating.The result output from the compiled program was in the form of a three-dimensional matrix providing the temperature at each node within the compact at various time steps.The data from the output files were plotted using MATLAB (ver.2008a).

Results and Discussion
Fig. 5 compares the effect of heating mode on the axial and radial thermal profile at various time steps for a high thermal conductivity cylindrical compact having radius of 0.65 cm.Fig.s 6a to 6c show the corresponding isothermal contours within cylindrical section half (Fig. 3b) for conventional heating, only microwave heating, and susceptor-assisted microwave heating, respectively.Tab.II summarizes the maximum temperature predicted from our simulation for various combinations of heating, compact size and its conductivity.
Tab. II.Effect of compact size and its conductivity on the maximum predicted temperature in conventional and microwave heating.The latter was simulated both with and without susceptor.
It is worth noticing that for high conductivity small-sized compacts, at the very initial stages (t=0.15s,0.30s), the temperature at the center of the sample is lower than that on the outside surface (top, bottom and side surface) (Fig. 5a).For conventional heating, the simulations were performed for ambient surface temperature of 1200K.It is worth noticing that for such small samples having high conductivity the temperature soon attains a steady state condition in just 1.8s (Fig. 6a).In comparison with conventional heating, during microwave heating the heat generation occurs within the sample per se.However, due to high conductivity of the sample, the temperature is uniform both across the radial and axial directions.However, even after 180s the maximum temperature attained within the compact is below 1200K.The only difference observed is that in the presence of the susceptor the maximum temperature attained is more than twice as much as that observed for similar timeframe when the compacts are heated only in microwaves.

Compact Size Thermal Conductivity
Heating Mode T max , K  The effect of heating mode on the simulated thermal profile and isothermal contours of a relatively low conductivity cylindrical compact to 0.65 cm radius is shown in Fig. s 7 and 8, respectively.As compared to the high conductivity counterparts, for similar dimensions, the low conductivity samples when heated conventionally (Figs.7a and 8a) show more pronounced thermal gradient for the first 10s.Consequently, it takes more time for the temperature to attain steady state condition.In contrast, for microwave heating (Fig. 7b) the samples temperature is much more uniform (Fig. 8b).The trend remains the same even in the presence of the susceptor.However, the maximum temperature attained is significantly higher in susceptor-assisted microwave heating (Figs.7c and 8c).Figs 9 and 10 show the effect of heating mode on the temporal thermal profile and contours for high conductivity compacts of 6.5 cm radius.In spite of the high conductivity, because of its relatively large size, conventional heating reveals steep temperature gradient in the beginning (Figs. 9a,10a).Consequently, to homogenize the temperature distribution one has to heat it for a longer duration.In actual systems, both these situations are detrimental to the properties [5].In case of microwave heating, the temperature is evenly distributed across the section (Figs.9b,10b).However, it does take about 180s for the temperature to reach 1200K.The same samples when subjected to microwaves in presence of SiC reveal very interesting trend (Fig. 9c).As evident from Table II and Fig. 10c the maximum temperature attained in the presence of the susceptor is much higher and even exceed that attained through conventional heating.This underscores the efficacy of susceptor-aided microwave heating.
Interestingly, initially, the susceptor couples more strongly than the sample.As a result the outside surface of the sample has a higher temperature.With time however the contribution of the same which too acts as a heat source becomes predominant.As a result at 180s, the temperature profile becomes inverse with the sample temperature at the center being higher than that attained outside.show the thermal profile and contours in a cylindrical compact of 6.5 cm having low conductivity heated in a conventional furnace.The thermal profile is remarkably different than those obtained for the high conductivity samples (Fig. 9a).As expected for low conductivity samples both the axial as well as radial thermal gradients are very steep.While the maximum temperature is 1200K, for the timeframe used in the simulation, the samples fail to attain a steady state condition.In contrast, for microwave heating condition (Figs.11b and 12b), the temperature attained is relatively uniform throughout the samples.For longer time, the compacts do tend to get an inverse thermal profile.Interestingly, in the presence of the susceptor the samples exhibit a mixed heating mode.For smaller time duration the temperature at the surface is higher than that in the inside region.For longer durations, the thermal profile becomes inverse in a more pronounced manner.Experimental results showing such a behavior has been reported by Binner and coworkers [28].From a design consideration, while it is recommended for medium-sized components -particularly, those having lower conductivity -to have susceptors, the time should be optimized so as to prevent the temperature runaways and inverse gradient.Figs 13 and 14 show the simulated temporal variation of temperatures within a high conductivity cylindrical samples having 65 cm radius along its radius and height.Note that despite the high conductivity, owing to the large sample size there is a thermal gradient during conventional heating.In fact, for all the heating modes, the behavior of this sample is akin to that shown by the medium sized sample (6.5 cm radius) having lower conductivity (Figs.11  and 12).As shown in Fig. 15a, for such large sized cylinders, the thermal gradients are much steeper when the samples have poor conductivity.In fact, for the same level of power setting and simulation time, the results prove that conventional heating is virtually ineffective with only the surface region being heated (Fig. 16a).This implies that for attaining the steady state condition the process time has to be significantly increased.In case of microwave heating however, the effective thermal mass of the sample which acts as the heat source is high.Hence, except the edges, the temperature remains homogenous throughout the sample (Figs.15b and 16b).In fact owing to the large sample size, the presence of susceptor does not significantly influence the temperature distribution within the compact (Figs.15c and 16c).This is also reflected in the maximum predicted temperature (Tab.II) which shows only 3.5% increase when microwave heating is done in conjunction with susceptor.3b) having radius of 0.65 cm, 6.5 cm and 65 cm, respectively.It is evident that irrespective of the conductivity, for small-sized compacts (0.65 cm) conventional heating is most effective.For medium sized compacts (6.5 cm) having higher thermal conductivity both conventional and microwave heating results in similar maximum temperature.In the presence of the susceptor though the temperature attained is higher.For low-conductivity medium-sized cylindrical samples (Fig. 17b) though conventional heating is hardly effective.In comparison, microwave (both with and without susceptor) result in significant enhancement of the sample temperature.The most significant difference is observed for large-sized compacts (Fig. 17c) wherein microwave are the most effective means to heat the samples in a short time span.The above results amply demonstrate the response of varying compact size and its conductivity to the various heating modes.It is envisaged that this will enable the design of an optimum heating methodology customized to the physical dimensions and thermal attributes of the specific system.

Conclusions
The present study investigates the effect of the sample dimension and its physical properties on the temporal temperature evolution in microwave heating.Using a 2-D finite difference approach the temperature distribution in cylindrical samples of varying radii (0.65, 6.5, 65 cm) and thermal conductivites was simulated for microwave as well as conventional heating.The influence of susceptor on the sample heating behaviour in microwaves was also modeled.In microwave heaing, it is possible to have an inverse thermal profile wherein the temperature at the core is higher than those achieved at the surface.The tendency to exhibit inverse thermal gradient increases with heating time and is higher for large sizes and lower conductivity samples.Use of susceptor aids is homogenization of temperature within the samples.For smaller-sized samples (0.65 cm) -irrespective of their thermal conductivitymicrowave heating is ineffective.In contrast, samples are uniformly heated to higher temperatures in conventional heating in a relatively short duration.To effectively heat such samples in microwaves susceptor is essential.For high conductivity medium-sized samples (6.5 cm) the heating behaviour is independent of the heating mode.The beneficial effect of microwave heating is realized for lower conductivity samples which further manifest itself in the presence of susceptor.For much larger sample size (65 cm) -irrespective of its thermal conductivity -microwave heating is a more practical option for attaining a uniform temperature within a short duration.For such large sample dimensions, susceptor does not influence the heating efficacy.Unlike microwave heating, conventional heating is restricted only to the sample surface.This implies a more prolonged heating period which may not always be the best option from the processing consideration.It is envisaged that the findings from the present study can be used as a guiding tool to customize optimum heating methodology that takes into consideration the thermal attributes and the physical dimensions of the sample.профили су симулирани за цилиндре са различитим презницима (0.65, 6.5 и 65 см) и за разлизита физичка својства.Резултати симулације указују на разлике и понашању узорака у микроталасима.Ефикасност микроталасног загревања зависи од величине честица и термалне проводности.Кључне речи: Микроталасно загревање, пренос топлоте, моделирање коначном разликом.

Fig. 1 .
Fig. 1.Schematic of the microwave setup used in the present study.

Fig. 2 .
Fig. 2. Flow diagram of the 2-D finite difference model used for cylindrical geometry.

Fig. 3 .
Fig. 3. (a) Discretization of the calculation domain of the cylindrical compact having the radius to height ratio of unity.(b) Highlighted area in represents the section for which axial and radial temperature profiles and isothermal contours were calculated.

Fig. 4 .
Fig. 4. Section of the cylindrical compact showing (a) the side surface and (b) top/bottom surface node.

Fig. 5 .
Fig. 5. Simulated axial and radial temperature distribution in a high thermal conductivity cylindrical compact having 0.65 cm radius under (a) conventional and MW heating without (b) and with (c) susceptor assistance.

Fig. 6 .
Fig. 6.Simulated temperature distribution showing isothermal contours as a function of time in a high conductivity cylindrical compact having 0.65 cm radius under conventional heating (a), MW heating without susceptor (b) and susceptor-assisted MW heating (c).

Fig. 7 .
Fig. 7. Simulated axial and radial temperature distribution in a low thermal conductivity cylindrical compact having 0.65 cm radius under (a) conventional and MW heating without (b) and with (c) susceptor assistance.

Fig. 8 .
Fig. 8. Simulated temperature distribution showing isothermal contours as a function of time in a low conductivity cylindrical compact having 0.65 cm radius under conventional heating (a), MW heating without susceptor (b) and susceptor-assisted MW heating (c).

Fig. 9 .
Fig. 9. Simulated axial and radial temperature distribution in a high thermal conductivity cylindrical compact having 6.5 cm radius under (a) conventional and MW heating without (b) and with (c) susceptor assistance.

Fig. 10 .
Fig. 10.Simulated temperature distribution showing isothermal contours as a function of time in a high conductivity cylindrical compact having 6.5 cm radius under conventional heating (a), MW heating without susceptor (b) and susceptor-assisted MW heating (c).

Fig. 11 .
Fig. 11.Simulated axial and radial temperature distribution in a low thermal conductivity cylindrical compact having 6.5 cm radius under (a) conventional and MW heating without (b) and with (c) susceptor assistance.

Fig.s 11a and 12a
Fig.s 11a and 12ashow the thermal profile and contours in a cylindrical compact of 6.5 cm having low conductivity heated in a conventional furnace.The thermal profile is remarkably different than those obtained for the high conductivity samples (Fig.9a).As expected for low conductivity samples both the axial as well as radial thermal gradients are very steep.While the maximum temperature is 1200K, for the timeframe used in the simulation, the samples fail to attain a steady state condition.In contrast, for microwave heating condition (Figs.11b and 12b), the temperature attained is relatively uniform throughout the samples.For longer time, the compacts do tend to get an inverse thermal profile.Interestingly, in the presence of the susceptor the samples exhibit a mixed heating mode.For smaller time duration the temperature at the surface is higher than that in the inside region.For longer durations, the thermal profile becomes inverse in a more pronounced manner.

Fig. 12 .
Fig. 12. Simulated temperature distribution showing isothermal contours as a function of time in a low conductivity cylindrical compact having 6.5 cm radius under conventional heating (a), MW heating without susceptor (b) and susceptor-assisted MW heating (c).

Fig. 13 .
Fig. 13.Simulated axial and radial temperature distribution in a high thermal conductivity cylindrical compact having 65 cm radius under (a) conventional and MW heating without (b) and with (c) susceptor assistance.

Fig. 14 .
Fig. 14.Simulated temperature distribution showing isothermal contours as a function of time in a high conductivity cylindrical compact having 65 cm radius under conventional heating (a), MW heating without susceptor (b) and susceptor-assisted MW heating (c).

Fig. 15 .
Fig. 15.Simulated axial and radial temperature distribution in a low thermal conductivity cylindrical compact having 65 cm radius under (a) conventional and MW heating without (b) and with (c) susceptor assistance.

Fig. 16 .
Fig. 16.Simulated temperature distribution showing isothermal contours as a function of time in a low conductivity cylindrical compact having 65 cm radius under conventional heating (a), MW heating without susceptor (b) and susceptor-assisted MW heating (c).

Fig. 17 .
Fig. 17.Effect of heating mode on the temporal thermal profile at the center of cylindrical compact (Fig. 3b) of radius (a) 0.65 cm, (b) 6.5 cm and (c) 65 cm having high (left) and low (right) thermal conductivity.

Fig.s 17a
Fig.s 17a to 17c compare the effect of heating mode and conductivity on the temperature variation at various time steps at the center of cylindrical compacts (Fig.3b) having radius of 0.65 cm, 6.5 cm and 65 cm, respectively.It is evident that irrespective of the conductivity, for small-sized compacts (0.65 cm) conventional heating is most effective.For medium sized compacts (6.5 cm) having higher thermal conductivity both conventional and microwave heating results in similar maximum temperature.In the presence of the susceptor