- Email: [email protected]
A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties
ARTICLE IN PRESS
JID: APM
[m3Gsc;November 18, 2016;21:36]
Applied Mathematical Modelling 0 0 0 (2016) 1–11
Contents lists available at ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties M.C. Navarro∗, J. Burgos Departamento de Matemáticas, Facultad de Ciencias y Tecnologías Químicas, Universidad de Castilla-La Mancha, Ciudad Real 13071, Spain
a r t i c l e
i n f o
Article history: Received 14 March 2016 Revised 5 October 2016 Accepted 27 October 2016 Available online xxx Keywords: Numerical modeling Microwave heating Spectral methods
a b s t r a c t In this paper, we develop a numerical model based on spectral methods for the simulation of heat transfer due to radial irradiation microwave applied to samples in cylindrical geometry. We solve the Maxwell’s equations and the resulting electric field distribution is incorporated as a source term in the heat transfer equation. The model includes the temperature dependence of the dielectric properties. The numerical model is validated with experimental temperature data from literature. © 2016 Elsevier Inc. All rights reserved.
1. Introduction Microwave heating is one of the most important methods for heating materials. It is used for different processes such as baking, thawing, heating, cooking and drying [1–5]. Microwaves offer several advantages: high speed startup, less process time, internal heating, energy efficiency, making no pollution, which have made them a high demand technology in industrial and household applications. Modeling of microwave heating involves coupling the models for microwave power absorption and temperature distribution inside the sample. Several simulation studies have modeled the heat generation due to microwaves by considering that the microwave power decreases exponentially as a function of penetration into the product. This simple approach, known as Lambert’s law, is valid only for large sample dimensions and high loss dielectric materials [6]. In other case, a complete solution of the unsteady Maxwell’s equations is required [7–9]. Simulation techniques have been extensively applied to model heat transfer due to microwaves. Among the computational methods available, Finite Element Method (FEM) and Finite Difference (FD) have been the most commonly used for solving microwave heating problems [10]. There are several studies on microwave heating of cylinders involving heating with radial irradiation using a finite element method [11–14]. The objective of this study is to provide a numerical method based on spectral methods for the simulation of heat transfer due to radial irradiation microwave in cylindrical geometry. Spectral methods are usually the best tool for solving partial differential equations to high accuracy on a simple domain if data defining the problem are smooth. They can achieve more accuracy than a finite difference or a finite element method. ∗
Corresponding author. Fax: +34 926295318. E-mail address: [email protected] (M.C. Navarro).
http://dx.doi.org/10.1016/j.apm.2016.10.062 0307-904X/© 2016 Elsevier Inc. All rights reserved.
Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM 2
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
Fig. 1. Physical setup. Electromagnetic radiation incident normal to the surface.
They also demand less computer memory than the alternatives. From the three most commonly used spectral schemes the collocation method has been chosen. This is a numerical method often used in thermoconvective problems [15–18] which has been demonstrated to be an efficient and useful tool for solving numerically the partial differential equations that model those problems. In this method, the approximation of a field w is given by its projection wN into a certain functional space of finite dimension (a space of polynomials). It is imposed that wN verifies the discretization of the operator we want to solve in certain points (collocation points), chosen depending on the problem and the selection of the basis for the space in which w is approximated. In the present work, we develop a numerical method based on a spectral collocation method to predict the temperature of a cylindrical sample heated radially by microwave.
2. Mathematical model formulation 2.1. Maxwell’s power dissipation The evaluation of the temperature distribution in any material subject to microwave irradiation depends on the knowledge of the electromagnetic field resulting from microwave power absorption. Several simulation studies model the heat generation due to microwaves by considering that the microwave power decreases exponentially as a function of penetration into the product. Such approach, known as Lambert’s law, can be obtained through a series of simplifications applied to Maxwell’s equations [19]. This is valid only for semi-infinite samples with dimensions much larger that the wave length. A more rigorous procedure to evaluate the electromagnetic field distribution consists on solving Maxwell’s equations:
∇ ×H=J+
∂D , ∂t
(1)
∇ ·D=q
(2)
∇ · B = 0,
(3)
where E is the electric field, B is the magnetic flux density, H is the magnetic field, J is the current density, D is the flux density and q the electric charge density. The constitutive relations are
J = σ E,
(4)
D = E,
(5)
B = μH ,
(6)
where σ is the electric conductivity, μ is the magnetic permeability and is the electric permittivity. With the help of the constitutive relations and assuming that the material magnetic permeability μ is approximated by its value in free space, considering electroneutral conditions of the material, and a one-dimensional analysis with the incident radiation assumed to be normal to the surface of the cylinder (Fig. 1), Maxwell’s equations yield:
d2 E 1 dE + + k21 E = 0, for 0 < r < R, r dr dr 2
(7)
Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
[m3Gsc;November 18, 2016;21:36] 3
where k21 = ω2 μ0 0 (κ + iκ ), being ω the angular velocity, μ0 the free space permeability, 0 the free space permittivity, κ the dielectric constant and κ the dielectric loss factor. After rescaling
r E and u = , R E0
r =
(8)
where E0 is the intensity of the incident field and u = v + iw, Eq. (7) is expressed as follows:
d2 v dr 2 d2 w dr 2
1 dv + + ψv − ξ w = 0, r dr
(9)
1 dw + + ψ w + ξ v = 0, r dr
(10)
where ψ = R2 ω2 μ0 0 κ , ξ = R2 ω2 μ0 0 κ . The boundary conditions at r = 0 are
dv dw = =0 dr dr
(11)
and boundary conditions at r = 1 are [20]:
dv + c1 v + c2 w = c3 , dr
(12)
dw + c1 w − c2 v = c4 , dr
(13)
where
c1 = Rα0
c2 =
c3 =
c4 =
J1 (α0 R )J0 (α0 R ) + Y1 (α0 R )Y0 (α0 R ) , J02 (α0 R ) + Y02 (α0 R ) 2
(14)
, π J02 (α0 R ) + Y02 (α0 R )
(15)
Y0 (α0 R ) , π J02 (α0 R ) + Y02 (α0 R )
(16)
J0 (α0 R ) , π J02 (α0 R ) + Y02 (α0 R )
(17)
−4
−4
being J and Y Bessel functions. Subscripts 0 and 1 refer to zero-order and first-order of the Bessel functions and α 0 is the free space wave number. Applying Poynting power theorem, the power dissipated per unit volume is given by:
Q (r ) =
1 ω0 κ E E ∗ , 2
(18)
where E∗ is the complex conjugate of E. 2.2. Heat equation and temperature distribution The heat transport within a solid sample due to volumetric absorption of microwave power Q is governed by the energy balance equation
ρC p
∂T = κ∇ 2 T + Q, ∂t
(19)
where T is the temperature, ρ Cp is the specific heat per unit volume and κ is the thermal conductivity. In cylindrical coordinates, equation (19) for a one-dimensional model becomes
ρC p ∂ T ∂ 2 T 1 ∂ T Q = + + , in 0 < r < R. κ ∂t ∂ r2 r ∂ r κ
(20)
The initial and boundary conditions considered are:
T = T0 at t = 0,
(21)
Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM 4
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
∂T = 0 at r = 0, t > 0, ∂r ∂T −κ = h(T − Ta ) at r = R, t > 0, ∂r
(22) (23)
where h is the heat transfer coefficient, Ta is the air temperature and T0 is the initial temperature of the sample. We have solved the problem in dimensionless form after rescaling r = r/R, t = κ t/(ρC p R2 ), T = κ (T − T0 )/(Q0 R2 ), Q = Q/Q0 , where Q0 is the initial power dissipated per unit volume. The dimensionless heat equation together with the initial and boundary conditions (with primes now omitted) become:
∂ T ∂ 2T 1 ∂ T = + + Q in 0 < r < 1, ∂t ∂ r2 r ∂ r
(24)
T = 0 at t = 0,
(25)
∂T = 0 at r = 0, t > 0, ∂r ∂ T hR h − = T+ (T0 − Ta ) at r = 1, t > 0. ∂r κ Q0 R
(26) (27)
3. Numerical implementation The schematic structure of the numerical method proposed consists on the following steps: 1. Compute κ (Tn ) and κ (Tn ) at temporal step n. 2. Solve Maxwell’s Eqs. (9)–(13) using the dielectric constant κ (Tn ) and the dielectric loss factor κ (Tn ) to obtain the electric field En at temporal step n. 3. Determine the heat generation Qn from Eq. (18). 4. Solve the heat equation (24) to obtain T n+1 , the temperature at step n + 1, using Qn and go back to 1. For the discretization of the heat equation we use an implicit Euler method at each temporal step:
T n+1 − T n = t
∂ 2 T n+1 1 ∂ T n+1 + + Q n in 0 < r < 1, r ∂r ∂ r2
T n+1 = 0 at t = 0,
∂ T n+1 = 0 at r = 0, t > 0, ∂r n+1 ∂T hR n+1 h − = T + (T0 − Ta ) at r = 1, t > 0. ∂r κ Q0 R
(28) (29) (30) (31)
We have used a dimensional time step t = 1 s in our computations. For the spatial discretization a Chebyshev collocation method is used. The unknown fields E and T at each temporal step are expanded in Chebyshev polynomials
E=
L−1
al Tl (r ),
(32)
bl Tl (r ),
(33)
l=0
T =
L−1 l=0
where Tn is the nth Chebyshev polynomial (see Appendix). Chebyshev polynomials are defined in the interval [−1, 1] and therefore, for computational convenience, the domain [0, 1] is transformed into [−1, 1]. This change of coordinates introduces scaling factors in equations and boundary conditions, which are not explicitly given here. There are L unknowns for E and T, which are determined by a collocation method. In particular, expansions (32) and (33) are substituted into the corresponding equation and boundary conditions and posed at the Gauss–Lobatto collocation points (see Appendix). At each temporal step linear systems of dimension L are found for coefficients al and bl, respectively. This can be easily solved with standard routines. Once the coefficients are found the fields are calculated from (32) and (33). The dimension used in the computations is L = 36. Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
JID: APM
ARTICLE IN PRESS M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
[m3Gsc;November 18, 2016;21:36] 5
Table 1 Thermal properties and input parameters for computations with agal gel. Parameter
Value
Density, ρ (kg/m ) Conductivity, k (W/(m °C)) Specific heat, Cp (J/(kg °C)) Heat transfer coefficient, h (W/m2 °C) Initial temperature, T0 (°C) Surrounding temperature, Ta (°C) Radius of cylinder, R (m) Dielectric constant, κ Dielectric loss, κ 3
10 0 0 0.5 4184 20 25 25 0.03 73.6 11.5
Fig. 2. (a) Temperature profiles predicted by the numerical model in present work (solid line) and experimental results obtained in Ref. [21] (symbols) after 15 s of microwave treatment; (b) contour for the predicted temperature.
4. Results and discussion 4.1. Validation of the mathematical model The numerical model developed in the present work has been validated with experimental temperature data from literature for microwave heating of foods and analogues. 4.1.1. Microwave heating of an agar gel cylinder The model has been first validated using data from [21,22] for a cylinder of agar gel of radius 0.03 m after 15 s of microwave treatment at 2450 MHz and 250 W. The dielectric properties are considered temperature independent. Thermal and electromagnetic properties and input parameters for the simulation are summarized in Table 1. In Fig. 2(a), experimental data by [21] and model predictions from our numerical simulations are shown. As observed, our model predicts correctly the typical temperature profile observed: microwaves focalize in the center, leading to hot spots and to uneven radial distribution of temperatures. It is noticed that the computational results show values below the experimental data at r = 0.01 m and r = 0.02 m. Despite this fact, the model predicts temperatures with high accuracy, obtained percentage relative errors between 1.4% and 3.9%. Fig. 2(b) shows the contour for the predicted temperature in the horizontal plane. In Fig. 3, experimental temperatures are shown together with the profile predicted for the center of the cylinder as a function of time. Experimental data used are reported in Ref. [22]. Predicted temperatures show good agreement with the behavior observed experimentally as a function of time. The average percentage relative error between predicted values and experimental data is 7.1%. It is appreciated a difference in the rate of heating at the center that gives a noticeable difference between the experimental and computed values at t = 60 s. It is remarkable that this difference is also found with other numerical models, as the finite difference method developed in Ref. [3]. The reason for this could be in the inputs used for the numerical computations. 4.1.2. Microwave heating of a mashed potato cylinder The first product tested in the case of dielectric properties dependent on the temperature is a cylinder of mashed potato. Thermal properties and input parameters for computations are summarized in Table 2. Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM 6
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
50 measured predicted
T (ºC)
40 30 20 10 0
0
20
40
60
80
100
120
140
160
180
t (s) Fig. 3. Comparison of the thermal history for a cylinder of agar gel at the center predicted by the numerical method in the present work (solid line) and experimental data (symbols) reported in Ref. [22] (T0 = 4◦ C, Ta = 60◦ C, h = 40 (W/m2 °C)). Table 2 Thermal properties and input parameters for computations with mashed potato. Parameter Density, ρ (kg/m3 ) Conductivity, k (W/(m °C)) Specific heat, Cp (J/(kg °C)) Heat transfer coefficient, h (W/m2 °C) Initial temperature, T0 (°C) Surrounding temperature, Ta (°C) Radius of cylinder, R (m) Dielectric constant, κ Dielectric loss, κ
Value 1056 0.55 3587 35 4 25 0.04 2.14 − 0.104 T + 66.822 3.09 − 0.0638 T + 17.615
a
b
Fig. 4. Temperature dependence of the dielectric constants for mashed potato at 2450 MHz. (a) Effect of temperature on the dielectric constant κ ; (b) effect of temperature on the dielectric loss factor κ .
The dependence of the dielectric constants with the temperature is linear as shown in Fig. 4. In Fig. 5(a), experimental data reported in Ref. [23] and model predictions (present work) are displayed for a cylinder of mashed potato of radius 0.04 m after 60 s of microwave treatment at 2450 MHz and 250 W. As observed, the numerical results agree with the experimental temperature profile observed. Temperature decreases from the surface to the center. The percentage relative errors found are between 0.96 and 4.6%, so the model predicts temperature accurately. Fig. 5(b) shows the contour for the predicted temperature in the horizontal plane showing the intense heating of the surface. 4.1.3. Microwave heating of a cylinder of kamaboko The second product tested is kamaboko, a traditional seafood in Japan that contains various kinds of surimi-based processes. Temperature distribution of kamaboko is analyzed in Ref. [24] experimentally and numerically by using a finite element method. The sample, made of surimi paste with 3 % NaCl and 5 % potato starch, is cylindrical and irradiated in Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
7
Fig. 5. (a) Temperature profiles predicted by the numerical model in present work (solid line) and experimental results (symbols) obtained in Ref. [23] in a cylinder of mashed potato after 60 s of microwave treatment at the 250-W oven setting, ; (b) contour for the predicted temperature. Table 3 Thermal properties and input parameters for computations with surimi. Parameter Density, ρ (kg/m3 ) Conductivity, k (W/(m° C)) Specific heat, Cp (J/(kg °C)) Heat transfer coefficient, h (W/m2 °C) Initial temperature, T0 (°C) Surrounding temperature, Ta (°C) Radius of cylinder, κ Dielectric constant, κ Dielectric loss, R(m)
Value 1050 0.55 3600 50 25.8 22 0.027 −0.0052 T 3 + 0.3624 T 2 − 10.4229 T + 152.9286 −0.0056 T 3 + 0.3804 T 2 − 10.1052 T + 120.2142
the radial direction. Thermal properties and input parameters for computations are summarized in Table 3. The dielectric properties are temperature dependent. Now the dependence is more complex than linear. Experimental data for the dependence of the dielectric constant κ and the dielectric loss factor κ can be found in Ref. [25]. We have used interpolating polynomials for those data. Fig. 6 displays the profile of the interpolating polynomials used for the dielectric constant and the dielectric loss factor for surimi (3 % of NaCl and 5 % potato starch) at 2450 MHz. In Fig. 7(a), experimental data reported in Ref. [24] and the numerical results from our model are shown for a cylinder of surimi of radius 0.027 m after 180 s of microwave treatment at 300 W and 2450 MHz. For the sample analyzed the numerical results show that the high temperature domain concentrates at the edge of the cylinder, as experimental data confirm, whereas in the central part the temperature is very low. The model predicts temperatures accurately, with percentage relative errors between 1% and 5%. The intense heating of the surface can be clearly seen in Fig. 7(b) that shows the contour for the predicted temperature in the horizontal plane. 4.2. Convergence of the numerical method To carry out a test on the convergence of the numerical method we have compared the solutions for different time steps and spatial expansions. 4.2.1. Variation of the time step t We have studied the influence of the time step t used on the numerical results obtained. Fig. 8 shows the temperature profiles for different dimensional time steps t and L = 36 for the spatial expansion. Fig. 8(a) displays the case of agar gel after 15 s of microwave heating. As observed, results are not significantly sensitive to the time step used, even when large values of t are considered. Relative differences between solutions with t = 1 and t = 0.1 are of order 10−3 . When the dielectric properties are temperature dependent, as for the examples of mashed potato and surimi, results are similar. Fig. 8(b) shows the temperature profiles obtained for different time steps t in the case of a cylinder of mashed potato after 60 s of microwave heating. Again, little dependence on the time step chosen is appreciated. Fig. 8(c) reports analogous results in the surimi case, for which the dependence of the dielectric properties on the temperature is more complex. Relative differences between solutions with t = 1 and t = 0.1 are of order 10−3 in both cases. Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM 8
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
a
b
Fig. 6. Temperature dependence of the dielectric properties of surimi at 2450 MHz and 3 % of NaCl and 5 % potato starch. (a) Effect of temperature on the dielectric constant κ ; (b) effect of temperature on the loss factor κ .
a
b
Fig. 7. (a) Temperature profiles predicted by the numerical model in present work (solid line) and experimental results (symbols) obtained in Ref. [24] for a cylinder of surimi of radius 0.027 m at 2450 MHz and 3 % of NaCl and 5 % potato starch after 180 s of microwave treatment at 300 W; (b) contour for the predicted temperature.
These results are remarkable as they show that the numerical spectral method developed is not significantly sensitive to the time step chosen, which permits to select large time steps, specially if the microwave heating time is larger, in order to reduce the computation time. 4.2.2. Variation of the spatial expansion L We have also studied the dependence of the solution on the spatial expansion. Fig. 9 shows the temperature profiles for different spatial expansions L = 16, L = 26 and L = 36, and a dimensional time step t = 1 for the three examples studied: agar gel, mashed potato and surimi. As observed in all cases, there is not a significant influence of the spatial expansion used on the results obtained. It is a characteristic of the spectral methods. With only 16 grid points the temperature profile is almost as accurate as that obtained with 36 grid points. Errors begin to be significant for L ≤ 10. The results reinforce the power of the numerical method developed in the present work. 4.2.3. Computational cost We have calculated the computational effort required to perform the simulations depending on the time step t and spatial expansion L. Table 4 shows the computational cost for the case of the mashed potato sample (60 s of microwave heating). The computational cost for the case of agar gel (15 s of microwave heating) and for the case of surimi (180 s of microwave heating) are essentially the values reported in Table 4 divided by four and multiplied by three, respectively. A 4GHz Intel Core i7 processor, with a RAM of 32 Gb 1867MHz has been used in our computations. As shown in Table 4, Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
a
b
9
c
Fig. 8. Predicted temperature profiles for different time steps t and L = 36 for the spatial expansion. (a) Agar gel after 15 s of microwave heating; (b) mashed potato after 60 s of microwave heating; (c) surimi after 180 s of microwave heating.
a
b
c
Fig. 9. Predicted temperature profiles for different spatial expansions L and a time step t = 1 and L = 36. (a) Agar gel after 15 s of microwave heating; (b) mashed potato after 60 s of microwave heating; (c) surimi after 180 s of microwave heating. Table 4 Computational cost (s) for the numerical simulations depending on the time step t and the spatial expansion L for the case of mashed potato after 60 s of microwave heating.
L = 12 L = 16 L = 26 L = 36
t = 0.1
t = 0.5
t = 1
t = 5
6.62 13.10 52.65 184.95
1.33 2.65 10.59 36.99
0.66 1.32 5.30 18.51
0.13 0.27 1.06 3.69
the computational cost grows exponentially as the time step is reduced and the spatial expansion increased. Therefore it is remarkable the ability of the method to get accurate solutions with a relatively large time step and small spatial expansion. t = 1 and L = 16, for example, give accurate computations for the case of mashed potato in less than 1.5 s. 5. Conclusions In this work, we have developed a numerical method based on a spectral collocation method to predict the temperature of a cylindrical sample heated radially by microwave. We solve the Maxwell’s equations and the resulting electric field distribution is incorporated as a source term in the heat transfer equation. The model includes the temperature dependence of the dielectric properties. The numerical model is validated with experimental temperature data from literature. The predicted results are in agreement with the experiments. The use of spectra methods permit to obtain accurate results with small spatial expansions and large time steps, which reinforces the efficiency of the numerical method developed. Acknowledgments This work was partially supported by the Research Grant MTM2015-68818-R MINECO (Spanish Government), which includes RDEF funds. Appendix. Chebyshev polynomials and collocation points The Chebyshev polynomials {Tk (x), k = 0, 1, . . .} are the eigenfunctions of the singular Sturm–Liouville problem
d dx
1 − x2
1/2
Tk (x ) + k2 1 − x2
−1/2
Tk (x ) = 0.
Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM 10
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
For any k, Tk (x) is even if k is even, and odd if k is odd. If Tk is such that Tk (1 ) = 1, then
Tk (x ) = cos (kθ )
θ = arccos x.
The trigonometric relation cos(k + 1 )θ + cos(k − 1 )θ = 2 cos θ cos kθ gives the recurrence relation
Tk (x ) = 2xTk−1 (x ) − Tk−2 (x ), k = 2, 3, . . . , with T0 (x ) = 1 and T1 (x ) = x. These polynomials are orthogonal with respect to the weighted scalar product
( f, g)ω =
1
−1
f gωdx,
− 1
where ω = 1 − x2 2 . Turning now to relations of interest for discrete Chebyshev series, explicit formulas for the quadrature points and weights are: 1. Chebyshev–Gauss
x j = cos
(2 j + 1 )π 2n + 2
,
ρj =
π n+1
, j = 0, . . . , n.
2. Chebyshev–Gauss–Radau
2π j x j = cos , j = 0, . . . , n, 2n + 1
ρj =
3. Chebyshev–Gauss–Lobatto
x j = cos
πj n
, j = 0, . . . , n,
ρj =
π
2n+1 2π 2n+2
si
j=0
si
1 ≤ j ≤ n.
π
2n
si
j = 0, n
π
si
1 ≤ j < n,
n
where xj are the points and ρ j the discrete weights. The most commonly used points are the Chebyshev–Gauss–Lobatto. If we define the discrete product (u, v)n as
(u, v )n =
n
v x j ρ j,
u xj
j=0
the following quadrature formula
(u, v )n = (u, v )ω is fulfilled exactly for P2n−δ with δ = −1, 0, 1 for the Chebyshev–Gauss, Chebyshev–Gauss–Radau y Chebyshev–Gauss– Lobattopoints,respectively, where P2n−δ is the space of polynomials of degree 2n − δ . References [1] A.K. Datta, R.C. Ananthswaran, Handbook of Microwave Technology for Food Applications, Marcel Dekker Inc., USA, 2001. [2] E. Turabi, G. Sumnu, S. Sahin, Optimization of baking of rice cakes in infrared-microwave combination oven by response surface methodology, Food Bioprocess Technol. 1 (2008) 64–73. [3] L.A. Campañone, N.E. Zaritzky, Mathematical analysis of microwave heating, J. Food Eng. 69 (2005) 359–368. [4] L.A. Campañone, N.E. Zaritzky, Mathematical modeling and simulation of microwave thawing of large solid foods under different operating conditions, J. Food Eng. 69 (2009) 359–368. [5] W. Cha-um, P. Rattanadecho, W. Padkee, Experimental and numerical analysis of microwave heating of water and oil using a rectangular wave guide: influence of sample sizes, positions, and microwave power, Food Bioprocess Technol. 4 (2011) 544–558. [6] Q. Zhang, T.H. Jackson, A. Ungan, Numerical modeling of microwave induced natural convection, Int. J. Heat Mass Transf. 43 (20 0 0) 2141–2154. [7] D.C. Dibben, A.C. Metaxas, Experimental verification of a finite element solution to heating problems in a multi-mode cavity, J. Microw. Power Electromagn. En. 29 (4) (1994) 242–251. [8] P. Verboten, A.K. Datta, N.T. Anh, N. Sheerlinck, B.M. Nicolai, Computation of airflow effects on heat and mass transfer in a microwave oven, J. Food Eng. 59 (2003) 181–190. [9] P. Rattanadecho, The simulation of microwave heating of wood using a rectangular wave guide: influence of frequency and sample size, Chem. Eng. Sci. 61 (2006) 4798–4811. [10] V.V. Yakovlev, Examination of contemporary electromagnetic software capable of modeling problems of microwave heating, in: M. Willert-Porada (Ed.), Advances in Microwaveand Radio Frequency Processing, Springer Verlag, Berlin Heidelberg, 2006, pp. 178–190. [11] B.J. Pangrle, K.G. Ayappa, H.T. Davis, E.A. Davis, J. Gordon, Microwave thawing of cylinders, A. I. Ch. E. J. 37 (12) (1991) 1789–1800. [12] M.E.C. Oliveira, A.S. Franca, Microwave heating of foodstuffs, J. Food Eng. 53 (4) (2002) 347–359. [13] M. Kostoglou, T.D. Karapantsios, Approximate computation of heat sources in axisymmetric microwave heating, A. I. Ch. E. J. 52 (1) (2006) 408–413. [14] T. Basak, Role of lateral and radial irradiations on efficient microwave processing of food cylinders, Chem. Eng. Sci. 62 (2007) 3185–3196. [15] M.C. Navarro, M.A. M., H. Herrero, Instabilities in buoyant flows under localized heating, Chaos 17 (2007) 023105. [16] M.C. Navarro, H. Herrero, Vortex generation by a convective instability in a cylindrical annulus non-homogeneously heated, Phys. D 40 (2011) 1181–1188. [17] M.C. Navarro, D. Castaño, H. Herrero, Thermoconvective instabilities to explain the main characteristics of a dust devil-like vortex, Phys. D 308 (2015) 109–115.
Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
JID: APM
ARTICLE IN PRESS M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
[m3Gsc;November 18, 2016;21:36] 11
[18] D. Castaño, M.C. Navarro, H. Herrero, Secondary whirls in thermoconvective vortices developed in a cylindrical annulus locally heated form below, Commun. Nonlin. Sci. Numer. Simul. 28 (2015) 201–209. [19] S.A. Barringer, E.A. Davis, J. Gordon, K.G. Ayappa, H.T. Davis, Microwave heating temperature profiles for thin slabs compared to Maxwell and Lambert law predictions, J. Food Sci. 60 (5) (1995) 1137–1142. [20] C.A. Balanis, Advanced Engineering Electromagnetics, Wiley, New York, 1989. [21] G.W. Padua, Microwave heating of agar gels containing sucrose, J. Food Sci. 58 (60) (1993) 1426–1428. [22] H.W. Yang, S. Gunasekaran, Temperature profiles in a cylindrical model food during pulsed microwave heating, J. Food Sci. 66 (7) (2002) 998–1004. [23] S. Gunasekaran, H.W. Yang, Optimization of pulsed microwave heating, J. Food Eng. 78 (2007) 1457–1462. [24] W. Mao, M. Watanabe, N. Sakai, Analysis of temperature distributions in kamaboko during microwave heating, J. Food Eng. 71 (2005) 187–192. [25] W. Mao, M. Watanabe, N. Sakai, Dielectric properties of surimi at 915 MHz and 2450 MHz as affected by temperature, salt and starch, Fish. Sci. 69 (2003) 1042–1047.
Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
JID: APM
[m3Gsc;November 18, 2016;21:36]
Applied Mathematical Modelling 0 0 0 (2016) 1–11
Contents lists available at ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties M.C. Navarro∗, J. Burgos Departamento de Matemáticas, Facultad de Ciencias y Tecnologías Químicas, Universidad de Castilla-La Mancha, Ciudad Real 13071, Spain
a r t i c l e
i n f o
Article history: Received 14 March 2016 Revised 5 October 2016 Accepted 27 October 2016 Available online xxx Keywords: Numerical modeling Microwave heating Spectral methods
a b s t r a c t In this paper, we develop a numerical model based on spectral methods for the simulation of heat transfer due to radial irradiation microwave applied to samples in cylindrical geometry. We solve the Maxwell’s equations and the resulting electric field distribution is incorporated as a source term in the heat transfer equation. The model includes the temperature dependence of the dielectric properties. The numerical model is validated with experimental temperature data from literature. © 2016 Elsevier Inc. All rights reserved.
1. Introduction Microwave heating is one of the most important methods for heating materials. It is used for different processes such as baking, thawing, heating, cooking and drying [1–5]. Microwaves offer several advantages: high speed startup, less process time, internal heating, energy efficiency, making no pollution, which have made them a high demand technology in industrial and household applications. Modeling of microwave heating involves coupling the models for microwave power absorption and temperature distribution inside the sample. Several simulation studies have modeled the heat generation due to microwaves by considering that the microwave power decreases exponentially as a function of penetration into the product. This simple approach, known as Lambert’s law, is valid only for large sample dimensions and high loss dielectric materials [6]. In other case, a complete solution of the unsteady Maxwell’s equations is required [7–9]. Simulation techniques have been extensively applied to model heat transfer due to microwaves. Among the computational methods available, Finite Element Method (FEM) and Finite Difference (FD) have been the most commonly used for solving microwave heating problems [10]. There are several studies on microwave heating of cylinders involving heating with radial irradiation using a finite element method [11–14]. The objective of this study is to provide a numerical method based on spectral methods for the simulation of heat transfer due to radial irradiation microwave in cylindrical geometry. Spectral methods are usually the best tool for solving partial differential equations to high accuracy on a simple domain if data defining the problem are smooth. They can achieve more accuracy than a finite difference or a finite element method. ∗
Corresponding author. Fax: +34 926295318. E-mail address: [email protected] (M.C. Navarro).
http://dx.doi.org/10.1016/j.apm.2016.10.062 0307-904X/© 2016 Elsevier Inc. All rights reserved.
Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM 2
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
Fig. 1. Physical setup. Electromagnetic radiation incident normal to the surface.
They also demand less computer memory than the alternatives. From the three most commonly used spectral schemes the collocation method has been chosen. This is a numerical method often used in thermoconvective problems [15–18] which has been demonstrated to be an efficient and useful tool for solving numerically the partial differential equations that model those problems. In this method, the approximation of a field w is given by its projection wN into a certain functional space of finite dimension (a space of polynomials). It is imposed that wN verifies the discretization of the operator we want to solve in certain points (collocation points), chosen depending on the problem and the selection of the basis for the space in which w is approximated. In the present work, we develop a numerical method based on a spectral collocation method to predict the temperature of a cylindrical sample heated radially by microwave.
2. Mathematical model formulation 2.1. Maxwell’s power dissipation The evaluation of the temperature distribution in any material subject to microwave irradiation depends on the knowledge of the electromagnetic field resulting from microwave power absorption. Several simulation studies model the heat generation due to microwaves by considering that the microwave power decreases exponentially as a function of penetration into the product. Such approach, known as Lambert’s law, can be obtained through a series of simplifications applied to Maxwell’s equations [19]. This is valid only for semi-infinite samples with dimensions much larger that the wave length. A more rigorous procedure to evaluate the electromagnetic field distribution consists on solving Maxwell’s equations:
∇ ×H=J+
∂D , ∂t
(1)
∇ ·D=q
(2)
∇ · B = 0,
(3)
where E is the electric field, B is the magnetic flux density, H is the magnetic field, J is the current density, D is the flux density and q the electric charge density. The constitutive relations are
J = σ E,
(4)
D = E,
(5)
B = μH ,
(6)
where σ is the electric conductivity, μ is the magnetic permeability and is the electric permittivity. With the help of the constitutive relations and assuming that the material magnetic permeability μ is approximated by its value in free space, considering electroneutral conditions of the material, and a one-dimensional analysis with the incident radiation assumed to be normal to the surface of the cylinder (Fig. 1), Maxwell’s equations yield:
d2 E 1 dE + + k21 E = 0, for 0 < r < R, r dr dr 2
(7)
Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
[m3Gsc;November 18, 2016;21:36] 3
where k21 = ω2 μ0 0 (κ + iκ ), being ω the angular velocity, μ0 the free space permeability, 0 the free space permittivity, κ the dielectric constant and κ the dielectric loss factor. After rescaling
r E and u = , R E0
r =
(8)
where E0 is the intensity of the incident field and u = v + iw, Eq. (7) is expressed as follows:
d2 v dr 2 d2 w dr 2
1 dv + + ψv − ξ w = 0, r dr
(9)
1 dw + + ψ w + ξ v = 0, r dr
(10)
where ψ = R2 ω2 μ0 0 κ , ξ = R2 ω2 μ0 0 κ . The boundary conditions at r = 0 are
dv dw = =0 dr dr
(11)
and boundary conditions at r = 1 are [20]:
dv + c1 v + c2 w = c3 , dr
(12)
dw + c1 w − c2 v = c4 , dr
(13)
where
c1 = Rα0
c2 =
c3 =
c4 =
J1 (α0 R )J0 (α0 R ) + Y1 (α0 R )Y0 (α0 R ) , J02 (α0 R ) + Y02 (α0 R ) 2
(14)
, π J02 (α0 R ) + Y02 (α0 R )
(15)
Y0 (α0 R ) , π J02 (α0 R ) + Y02 (α0 R )
(16)
J0 (α0 R ) , π J02 (α0 R ) + Y02 (α0 R )
(17)
−4
−4
being J and Y Bessel functions. Subscripts 0 and 1 refer to zero-order and first-order of the Bessel functions and α 0 is the free space wave number. Applying Poynting power theorem, the power dissipated per unit volume is given by:
Q (r ) =
1 ω0 κ E E ∗ , 2
(18)
where E∗ is the complex conjugate of E. 2.2. Heat equation and temperature distribution The heat transport within a solid sample due to volumetric absorption of microwave power Q is governed by the energy balance equation
ρC p
∂T = κ∇ 2 T + Q, ∂t
(19)
where T is the temperature, ρ Cp is the specific heat per unit volume and κ is the thermal conductivity. In cylindrical coordinates, equation (19) for a one-dimensional model becomes
ρC p ∂ T ∂ 2 T 1 ∂ T Q = + + , in 0 < r < R. κ ∂t ∂ r2 r ∂ r κ
(20)
The initial and boundary conditions considered are:
T = T0 at t = 0,
(21)
Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM 4
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
∂T = 0 at r = 0, t > 0, ∂r ∂T −κ = h(T − Ta ) at r = R, t > 0, ∂r
(22) (23)
where h is the heat transfer coefficient, Ta is the air temperature and T0 is the initial temperature of the sample. We have solved the problem in dimensionless form after rescaling r = r/R, t = κ t/(ρC p R2 ), T = κ (T − T0 )/(Q0 R2 ), Q = Q/Q0 , where Q0 is the initial power dissipated per unit volume. The dimensionless heat equation together with the initial and boundary conditions (with primes now omitted) become:
∂ T ∂ 2T 1 ∂ T = + + Q in 0 < r < 1, ∂t ∂ r2 r ∂ r
(24)
T = 0 at t = 0,
(25)
∂T = 0 at r = 0, t > 0, ∂r ∂ T hR h − = T+ (T0 − Ta ) at r = 1, t > 0. ∂r κ Q0 R
(26) (27)
3. Numerical implementation The schematic structure of the numerical method proposed consists on the following steps: 1. Compute κ (Tn ) and κ (Tn ) at temporal step n. 2. Solve Maxwell’s Eqs. (9)–(13) using the dielectric constant κ (Tn ) and the dielectric loss factor κ (Tn ) to obtain the electric field En at temporal step n. 3. Determine the heat generation Qn from Eq. (18). 4. Solve the heat equation (24) to obtain T n+1 , the temperature at step n + 1, using Qn and go back to 1. For the discretization of the heat equation we use an implicit Euler method at each temporal step:
T n+1 − T n = t
∂ 2 T n+1 1 ∂ T n+1 + + Q n in 0 < r < 1, r ∂r ∂ r2
T n+1 = 0 at t = 0,
∂ T n+1 = 0 at r = 0, t > 0, ∂r n+1 ∂T hR n+1 h − = T + (T0 − Ta ) at r = 1, t > 0. ∂r κ Q0 R
(28) (29) (30) (31)
We have used a dimensional time step t = 1 s in our computations. For the spatial discretization a Chebyshev collocation method is used. The unknown fields E and T at each temporal step are expanded in Chebyshev polynomials
E=
L−1
al Tl (r ),
(32)
bl Tl (r ),
(33)
l=0
T =
L−1 l=0
where Tn is the nth Chebyshev polynomial (see Appendix). Chebyshev polynomials are defined in the interval [−1, 1] and therefore, for computational convenience, the domain [0, 1] is transformed into [−1, 1]. This change of coordinates introduces scaling factors in equations and boundary conditions, which are not explicitly given here. There are L unknowns for E and T, which are determined by a collocation method. In particular, expansions (32) and (33) are substituted into the corresponding equation and boundary conditions and posed at the Gauss–Lobatto collocation points (see Appendix). At each temporal step linear systems of dimension L are found for coefficients al and bl, respectively. This can be easily solved with standard routines. Once the coefficients are found the fields are calculated from (32) and (33). The dimension used in the computations is L = 36. Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
JID: APM
ARTICLE IN PRESS M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
[m3Gsc;November 18, 2016;21:36] 5
Table 1 Thermal properties and input parameters for computations with agal gel. Parameter
Value
Density, ρ (kg/m ) Conductivity, k (W/(m °C)) Specific heat, Cp (J/(kg °C)) Heat transfer coefficient, h (W/m2 °C) Initial temperature, T0 (°C) Surrounding temperature, Ta (°C) Radius of cylinder, R (m) Dielectric constant, κ Dielectric loss, κ 3
10 0 0 0.5 4184 20 25 25 0.03 73.6 11.5
Fig. 2. (a) Temperature profiles predicted by the numerical model in present work (solid line) and experimental results obtained in Ref. [21] (symbols) after 15 s of microwave treatment; (b) contour for the predicted temperature.
4. Results and discussion 4.1. Validation of the mathematical model The numerical model developed in the present work has been validated with experimental temperature data from literature for microwave heating of foods and analogues. 4.1.1. Microwave heating of an agar gel cylinder The model has been first validated using data from [21,22] for a cylinder of agar gel of radius 0.03 m after 15 s of microwave treatment at 2450 MHz and 250 W. The dielectric properties are considered temperature independent. Thermal and electromagnetic properties and input parameters for the simulation are summarized in Table 1. In Fig. 2(a), experimental data by [21] and model predictions from our numerical simulations are shown. As observed, our model predicts correctly the typical temperature profile observed: microwaves focalize in the center, leading to hot spots and to uneven radial distribution of temperatures. It is noticed that the computational results show values below the experimental data at r = 0.01 m and r = 0.02 m. Despite this fact, the model predicts temperatures with high accuracy, obtained percentage relative errors between 1.4% and 3.9%. Fig. 2(b) shows the contour for the predicted temperature in the horizontal plane. In Fig. 3, experimental temperatures are shown together with the profile predicted for the center of the cylinder as a function of time. Experimental data used are reported in Ref. [22]. Predicted temperatures show good agreement with the behavior observed experimentally as a function of time. The average percentage relative error between predicted values and experimental data is 7.1%. It is appreciated a difference in the rate of heating at the center that gives a noticeable difference between the experimental and computed values at t = 60 s. It is remarkable that this difference is also found with other numerical models, as the finite difference method developed in Ref. [3]. The reason for this could be in the inputs used for the numerical computations. 4.1.2. Microwave heating of a mashed potato cylinder The first product tested in the case of dielectric properties dependent on the temperature is a cylinder of mashed potato. Thermal properties and input parameters for computations are summarized in Table 2. Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM 6
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
50 measured predicted
T (ºC)
40 30 20 10 0
0
20
40
60
80
100
120
140
160
180
t (s) Fig. 3. Comparison of the thermal history for a cylinder of agar gel at the center predicted by the numerical method in the present work (solid line) and experimental data (symbols) reported in Ref. [22] (T0 = 4◦ C, Ta = 60◦ C, h = 40 (W/m2 °C)). Table 2 Thermal properties and input parameters for computations with mashed potato. Parameter Density, ρ (kg/m3 ) Conductivity, k (W/(m °C)) Specific heat, Cp (J/(kg °C)) Heat transfer coefficient, h (W/m2 °C) Initial temperature, T0 (°C) Surrounding temperature, Ta (°C) Radius of cylinder, R (m) Dielectric constant, κ Dielectric loss, κ
Value 1056 0.55 3587 35 4 25 0.04 2.14 − 0.104 T + 66.822 3.09 − 0.0638 T + 17.615
a
b
Fig. 4. Temperature dependence of the dielectric constants for mashed potato at 2450 MHz. (a) Effect of temperature on the dielectric constant κ ; (b) effect of temperature on the dielectric loss factor κ .
The dependence of the dielectric constants with the temperature is linear as shown in Fig. 4. In Fig. 5(a), experimental data reported in Ref. [23] and model predictions (present work) are displayed for a cylinder of mashed potato of radius 0.04 m after 60 s of microwave treatment at 2450 MHz and 250 W. As observed, the numerical results agree with the experimental temperature profile observed. Temperature decreases from the surface to the center. The percentage relative errors found are between 0.96 and 4.6%, so the model predicts temperature accurately. Fig. 5(b) shows the contour for the predicted temperature in the horizontal plane showing the intense heating of the surface. 4.1.3. Microwave heating of a cylinder of kamaboko The second product tested is kamaboko, a traditional seafood in Japan that contains various kinds of surimi-based processes. Temperature distribution of kamaboko is analyzed in Ref. [24] experimentally and numerically by using a finite element method. The sample, made of surimi paste with 3 % NaCl and 5 % potato starch, is cylindrical and irradiated in Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
7
Fig. 5. (a) Temperature profiles predicted by the numerical model in present work (solid line) and experimental results (symbols) obtained in Ref. [23] in a cylinder of mashed potato after 60 s of microwave treatment at the 250-W oven setting, ; (b) contour for the predicted temperature. Table 3 Thermal properties and input parameters for computations with surimi. Parameter Density, ρ (kg/m3 ) Conductivity, k (W/(m° C)) Specific heat, Cp (J/(kg °C)) Heat transfer coefficient, h (W/m2 °C) Initial temperature, T0 (°C) Surrounding temperature, Ta (°C) Radius of cylinder, κ Dielectric constant, κ Dielectric loss, R(m)
Value 1050 0.55 3600 50 25.8 22 0.027 −0.0052 T 3 + 0.3624 T 2 − 10.4229 T + 152.9286 −0.0056 T 3 + 0.3804 T 2 − 10.1052 T + 120.2142
the radial direction. Thermal properties and input parameters for computations are summarized in Table 3. The dielectric properties are temperature dependent. Now the dependence is more complex than linear. Experimental data for the dependence of the dielectric constant κ and the dielectric loss factor κ can be found in Ref. [25]. We have used interpolating polynomials for those data. Fig. 6 displays the profile of the interpolating polynomials used for the dielectric constant and the dielectric loss factor for surimi (3 % of NaCl and 5 % potato starch) at 2450 MHz. In Fig. 7(a), experimental data reported in Ref. [24] and the numerical results from our model are shown for a cylinder of surimi of radius 0.027 m after 180 s of microwave treatment at 300 W and 2450 MHz. For the sample analyzed the numerical results show that the high temperature domain concentrates at the edge of the cylinder, as experimental data confirm, whereas in the central part the temperature is very low. The model predicts temperatures accurately, with percentage relative errors between 1% and 5%. The intense heating of the surface can be clearly seen in Fig. 7(b) that shows the contour for the predicted temperature in the horizontal plane. 4.2. Convergence of the numerical method To carry out a test on the convergence of the numerical method we have compared the solutions for different time steps and spatial expansions. 4.2.1. Variation of the time step t We have studied the influence of the time step t used on the numerical results obtained. Fig. 8 shows the temperature profiles for different dimensional time steps t and L = 36 for the spatial expansion. Fig. 8(a) displays the case of agar gel after 15 s of microwave heating. As observed, results are not significantly sensitive to the time step used, even when large values of t are considered. Relative differences between solutions with t = 1 and t = 0.1 are of order 10−3 . When the dielectric properties are temperature dependent, as for the examples of mashed potato and surimi, results are similar. Fig. 8(b) shows the temperature profiles obtained for different time steps t in the case of a cylinder of mashed potato after 60 s of microwave heating. Again, little dependence on the time step chosen is appreciated. Fig. 8(c) reports analogous results in the surimi case, for which the dependence of the dielectric properties on the temperature is more complex. Relative differences between solutions with t = 1 and t = 0.1 are of order 10−3 in both cases. Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM 8
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
a
b
Fig. 6. Temperature dependence of the dielectric properties of surimi at 2450 MHz and 3 % of NaCl and 5 % potato starch. (a) Effect of temperature on the dielectric constant κ ; (b) effect of temperature on the loss factor κ .
a
b
Fig. 7. (a) Temperature profiles predicted by the numerical model in present work (solid line) and experimental results (symbols) obtained in Ref. [24] for a cylinder of surimi of radius 0.027 m at 2450 MHz and 3 % of NaCl and 5 % potato starch after 180 s of microwave treatment at 300 W; (b) contour for the predicted temperature.
These results are remarkable as they show that the numerical spectral method developed is not significantly sensitive to the time step chosen, which permits to select large time steps, specially if the microwave heating time is larger, in order to reduce the computation time. 4.2.2. Variation of the spatial expansion L We have also studied the dependence of the solution on the spatial expansion. Fig. 9 shows the temperature profiles for different spatial expansions L = 16, L = 26 and L = 36, and a dimensional time step t = 1 for the three examples studied: agar gel, mashed potato and surimi. As observed in all cases, there is not a significant influence of the spatial expansion used on the results obtained. It is a characteristic of the spectral methods. With only 16 grid points the temperature profile is almost as accurate as that obtained with 36 grid points. Errors begin to be significant for L ≤ 10. The results reinforce the power of the numerical method developed in the present work. 4.2.3. Computational cost We have calculated the computational effort required to perform the simulations depending on the time step t and spatial expansion L. Table 4 shows the computational cost for the case of the mashed potato sample (60 s of microwave heating). The computational cost for the case of agar gel (15 s of microwave heating) and for the case of surimi (180 s of microwave heating) are essentially the values reported in Table 4 divided by four and multiplied by three, respectively. A 4GHz Intel Core i7 processor, with a RAM of 32 Gb 1867MHz has been used in our computations. As shown in Table 4, Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
a
b
9
c
Fig. 8. Predicted temperature profiles for different time steps t and L = 36 for the spatial expansion. (a) Agar gel after 15 s of microwave heating; (b) mashed potato after 60 s of microwave heating; (c) surimi after 180 s of microwave heating.
a
b
c
Fig. 9. Predicted temperature profiles for different spatial expansions L and a time step t = 1 and L = 36. (a) Agar gel after 15 s of microwave heating; (b) mashed potato after 60 s of microwave heating; (c) surimi after 180 s of microwave heating. Table 4 Computational cost (s) for the numerical simulations depending on the time step t and the spatial expansion L for the case of mashed potato after 60 s of microwave heating.
L = 12 L = 16 L = 26 L = 36
t = 0.1
t = 0.5
t = 1
t = 5
6.62 13.10 52.65 184.95
1.33 2.65 10.59 36.99
0.66 1.32 5.30 18.51
0.13 0.27 1.06 3.69
the computational cost grows exponentially as the time step is reduced and the spatial expansion increased. Therefore it is remarkable the ability of the method to get accurate solutions with a relatively large time step and small spatial expansion. t = 1 and L = 16, for example, give accurate computations for the case of mashed potato in less than 1.5 s. 5. Conclusions In this work, we have developed a numerical method based on a spectral collocation method to predict the temperature of a cylindrical sample heated radially by microwave. We solve the Maxwell’s equations and the resulting electric field distribution is incorporated as a source term in the heat transfer equation. The model includes the temperature dependence of the dielectric properties. The numerical model is validated with experimental temperature data from literature. The predicted results are in agreement with the experiments. The use of spectra methods permit to obtain accurate results with small spatial expansions and large time steps, which reinforces the efficiency of the numerical method developed. Acknowledgments This work was partially supported by the Research Grant MTM2015-68818-R MINECO (Spanish Government), which includes RDEF funds. Appendix. Chebyshev polynomials and collocation points The Chebyshev polynomials {Tk (x), k = 0, 1, . . .} are the eigenfunctions of the singular Sturm–Liouville problem
d dx
1 − x2
1/2
Tk (x ) + k2 1 − x2
−1/2
Tk (x ) = 0.
Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
ARTICLE IN PRESS
JID: APM 10
[m3Gsc;November 18, 2016;21:36]
M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
For any k, Tk (x) is even if k is even, and odd if k is odd. If Tk is such that Tk (1 ) = 1, then
Tk (x ) = cos (kθ )
θ = arccos x.
The trigonometric relation cos(k + 1 )θ + cos(k − 1 )θ = 2 cos θ cos kθ gives the recurrence relation
Tk (x ) = 2xTk−1 (x ) − Tk−2 (x ), k = 2, 3, . . . , with T0 (x ) = 1 and T1 (x ) = x. These polynomials are orthogonal with respect to the weighted scalar product
( f, g)ω =
1
−1
f gωdx,
− 1
where ω = 1 − x2 2 . Turning now to relations of interest for discrete Chebyshev series, explicit formulas for the quadrature points and weights are: 1. Chebyshev–Gauss
x j = cos
(2 j + 1 )π 2n + 2
,
ρj =
π n+1
, j = 0, . . . , n.
2. Chebyshev–Gauss–Radau
2π j x j = cos , j = 0, . . . , n, 2n + 1
ρj =
3. Chebyshev–Gauss–Lobatto
x j = cos
πj n
, j = 0, . . . , n,
ρj =
π
2n+1 2π 2n+2
si
j=0
si
1 ≤ j ≤ n.
π
2n
si
j = 0, n
π
si
1 ≤ j < n,
n
where xj are the points and ρ j the discrete weights. The most commonly used points are the Chebyshev–Gauss–Lobatto. If we define the discrete product (u, v)n as
(u, v )n =
n
v x j ρ j,
u xj
j=0
the following quadrature formula
(u, v )n = (u, v )ω is fulfilled exactly for P2n−δ with δ = −1, 0, 1 for the Chebyshev–Gauss, Chebyshev–Gauss–Radau y Chebyshev–Gauss– Lobattopoints,respectively, where P2n−δ is the space of polynomials of degree 2n − δ . References [1] A.K. Datta, R.C. Ananthswaran, Handbook of Microwave Technology for Food Applications, Marcel Dekker Inc., USA, 2001. [2] E. Turabi, G. Sumnu, S. Sahin, Optimization of baking of rice cakes in infrared-microwave combination oven by response surface methodology, Food Bioprocess Technol. 1 (2008) 64–73. [3] L.A. Campañone, N.E. Zaritzky, Mathematical analysis of microwave heating, J. Food Eng. 69 (2005) 359–368. [4] L.A. Campañone, N.E. Zaritzky, Mathematical modeling and simulation of microwave thawing of large solid foods under different operating conditions, J. Food Eng. 69 (2009) 359–368. [5] W. Cha-um, P. Rattanadecho, W. Padkee, Experimental and numerical analysis of microwave heating of water and oil using a rectangular wave guide: influence of sample sizes, positions, and microwave power, Food Bioprocess Technol. 4 (2011) 544–558. [6] Q. Zhang, T.H. Jackson, A. Ungan, Numerical modeling of microwave induced natural convection, Int. J. Heat Mass Transf. 43 (20 0 0) 2141–2154. [7] D.C. Dibben, A.C. Metaxas, Experimental verification of a finite element solution to heating problems in a multi-mode cavity, J. Microw. Power Electromagn. En. 29 (4) (1994) 242–251. [8] P. Verboten, A.K. Datta, N.T. Anh, N. Sheerlinck, B.M. Nicolai, Computation of airflow effects on heat and mass transfer in a microwave oven, J. Food Eng. 59 (2003) 181–190. [9] P. Rattanadecho, The simulation of microwave heating of wood using a rectangular wave guide: influence of frequency and sample size, Chem. Eng. Sci. 61 (2006) 4798–4811. [10] V.V. Yakovlev, Examination of contemporary electromagnetic software capable of modeling problems of microwave heating, in: M. Willert-Porada (Ed.), Advances in Microwaveand Radio Frequency Processing, Springer Verlag, Berlin Heidelberg, 2006, pp. 178–190. [11] B.J. Pangrle, K.G. Ayappa, H.T. Davis, E.A. Davis, J. Gordon, Microwave thawing of cylinders, A. I. Ch. E. J. 37 (12) (1991) 1789–1800. [12] M.E.C. Oliveira, A.S. Franca, Microwave heating of foodstuffs, J. Food Eng. 53 (4) (2002) 347–359. [13] M. Kostoglou, T.D. Karapantsios, Approximate computation of heat sources in axisymmetric microwave heating, A. I. Ch. E. J. 52 (1) (2006) 408–413. [14] T. Basak, Role of lateral and radial irradiations on efficient microwave processing of food cylinders, Chem. Eng. Sci. 62 (2007) 3185–3196. [15] M.C. Navarro, M.A. M., H. Herrero, Instabilities in buoyant flows under localized heating, Chaos 17 (2007) 023105. [16] M.C. Navarro, H. Herrero, Vortex generation by a convective instability in a cylindrical annulus non-homogeneously heated, Phys. D 40 (2011) 1181–1188. [17] M.C. Navarro, D. Castaño, H. Herrero, Thermoconvective instabilities to explain the main characteristics of a dust devil-like vortex, Phys. D 308 (2015) 109–115.
Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062
JID: APM
ARTICLE IN PRESS M.C. Navarro, J. Burgos / Applied Mathematical Modelling 000 (2016) 1–11
[m3Gsc;November 18, 2016;21:36] 11
[18] D. Castaño, M.C. Navarro, H. Herrero, Secondary whirls in thermoconvective vortices developed in a cylindrical annulus locally heated form below, Commun. Nonlin. Sci. Numer. Simul. 28 (2015) 201–209. [19] S.A. Barringer, E.A. Davis, J. Gordon, K.G. Ayappa, H.T. Davis, Microwave heating temperature profiles for thin slabs compared to Maxwell and Lambert law predictions, J. Food Sci. 60 (5) (1995) 1137–1142. [20] C.A. Balanis, Advanced Engineering Electromagnetics, Wiley, New York, 1989. [21] G.W. Padua, Microwave heating of agar gels containing sucrose, J. Food Sci. 58 (60) (1993) 1426–1428. [22] H.W. Yang, S. Gunasekaran, Temperature profiles in a cylindrical model food during pulsed microwave heating, J. Food Sci. 66 (7) (2002) 998–1004. [23] S. Gunasekaran, H.W. Yang, Optimization of pulsed microwave heating, J. Food Eng. 78 (2007) 1457–1462. [24] W. Mao, M. Watanabe, N. Sakai, Analysis of temperature distributions in kamaboko during microwave heating, J. Food Eng. 71 (2005) 187–192. [25] W. Mao, M. Watanabe, N. Sakai, Dielectric properties of surimi at 915 MHz and 2450 MHz as affected by temperature, salt and starch, Fish. Sci. 69 (2003) 1042–1047.
Please cite this article as: M.C. Navarro, J. Burgos, A spectral method for numerical modeling of radial microwave heating in cylindrical samples with temperature dependent dielectric properties, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.10.062