Numerical Investigation of the Deformable Porous Media Treated by the Intermittent Microwave

: A 2D axi-symmetric theoretical model of dielectric porous media in intermittent microwave (IMW) thermal process was developed, and the electromagnetic energy, multiphase transport, phase change, large deformation, and glass transition were taken into consideration. From the simulation results, the mass was mainly carried by the liquid water, and the heat was mainly carried by liquid water and solid. The diffusion was the dominant mechanism of the mass transport during the whole process, whereas for the heat transport, the convection dominated the heat transport near the surface areas during the heating stage. The von Mises stress reached local maxima at different locations at different stages, and all were lower than the fracture stress. A material treated by a longer intermittent cycle length with the same pulse ratio (PR) tended to trigger the phenomena of overheat and fracture due to the more intense ﬂuctuation of moisture content, temperature, deformation, and von Mises stress. The model can be extended to simulate the intermittent radio frequency (IRF) process on the basis of which one can select a suitable energy source for a speciﬁc process.


Introduction
Microwave energy, classified as an electromagnetic energy, is one of the promising green resources, also having the advantages of high efficiency and rapid volumetric heating compared with the conventional heating methods. In order to prevent interfering the telecommunications, special frequency bands are reserved for industrial, scientific, and medical applications (ISM). The allowed microwave frequency bands are located at 433 MHz, 915 MHz and 2450 MHz [1]. The thermal effect of microwave has been applied in various processes of dielectric materials (e.g., drying [2,3], pasteurization [4], disinsection [5]). However, the non-uniform electric field is the major limitation of microwave applications. In order to overcome this drawback, the usage of the intermittent microwave source is one of the effective solutions [6]. Apart from microwave energy, a commonly used electromagnetic energy source is radio frequency, whose permitted frequency bands used in ISM applications are located at 13.56 MHz, 27.12 MHz, and 40.68 MHz [7]. The differences between these two electromagnetic energy sources are the structure of the equipment, the deep penetration, the heating efficiency, etc. [8], which are all caused by the frequency differences. Therefore, the two electromagnetic energies share many characters and can be studied and discussed together.
Most of the dielectric materials can be treated as porous media, such as food materials [9,10], minerals and fuels [11,12], clay and ceramics [13], and animal and plant tissues [14,15]. During the thermal treatments of the porous dielectric materials by the microwave energy, the mass and heat within the products tend to be transported under the gradients of concentration, temperature, and pressure [16][17][18], and sequentially, internal stresses, structure collapses, volumetric changes, or even crack and fracture would occur [11,19]. In order to reveal the mechanism and optimize the complex processes, researchers have carried out many studies. Experiments coupled with various measuring technologies are the basis of simulation studies. They provide the physical validation for the simulation studies, obtaining critical empirical property parameters. However, the setups of the experiments and the subsequent measuring instruments are costly in terms of time and money. Thus, as a supplement, physically based theoretical models provide the direct description of various thermal processes efficiently, which guide the optimization of the experiments in turn [20,21].
In order to reduce the complexity and increase the computability of the theoretical models, some simplifying assumptions are inevitable. The simplest one is that the effects of deformation are introduced by input parameters obtained experimentally or empirically, e.g., the shrinkage dependent diffusivity [22,23] or the shrinkage velocity [24,25]. Furthermore, the single mechanical deformation is considered in some theoretical models, e.g., the hyperelastic deformation [26][27][28] and viscoelastic deformation [29]. However, different samples during diverse processes perform not only the mechanical deformation but also the deformation due to the moisture or temperature changes. Namely, when the material deformation is small, the total deformation is the summation of the mechanical deformation (e.g., elastic deformation [30,31], elastoplastic deformation [32,33], viscoelastic deformation [34][35][36][37]) and the deformation due to the moisture and temperature changes [11,[38][39][40]. While for a large deformation, the total deformation is multiplicative, split into the mechanical deformation (e.g., hyperelastic deformation [41][42][43]) and the moisture loss-induced deformation.
During the calculation of the mechanical deformation, Young's modulus and Poisson's ratio are two critical parameters that directly influence the results of deformation and stress. From the previous studies, Young's modulus is used as a constant [31,[34][35][36]44] or a function of moisture content [30,32,33,42,43,45]. In most of the theoretical models, the materials are assumed to be solely rubber-like and undergo an ideal shrinkage during the processes [32,33,36,41,43,45]. However, the Poisson's ratio, a more realistic description, is considered as a variable, varying with moisture and temperature [46][47][48]. This is because some materials will turn from a soft and rubbery-like state to a hard and glass-like state when the moisture content decreases or the temperature increases [42,45], and vice versa [49,50].
The purpose of this research was to develop a comprehensive physically based theoretical model of the intermittent microwave (IMW) thermal process that takes the electromagnetic energy, multiphase transport, phase change, large deformation, and glass transition into consideration. The theoretical model is more flexible than the usage of experiments or the dynamic models. Thus, various industrial thermal processes of different dielectric porous media can be explained in depth and be optimized in an efficient and economical way. For example, the theoretical model could be used to cope with the parameter optimization of various microwave thermal processes (e.g., drying, heating, thawing, blanching) of different porous media (e.g., food, wood, ceramic, minerals). Meanwhile, the temperature, moisture content, and the deformation obtained by the model could be used as the evaluation indexes. Therefore, the cost in terms of time, manpower, and resources could be saved.
In this manuscript, the simulation results of moisture content, deformation, and temperature were validated by experimental data of the IMW thermal process. Moreover, the developed model was used to simulate the process in the reference [18]. Discussions about the mechanisms of the mass and heat transport were induced by introducing the mass and heat Peclet (Pe) numbers. Further discussions about the effects of the intermittent cycle length with the same pulse ratio (PR) were carried out. Finally, the model was extended to simulate the intermittent radio frequency (IRF) thermal process, and the results of the IMW thermal process and the IRF thermal process were compared and discussed.

Process Description and Assumptions
As shown in Figure 1a, a hygroscopic cylindrical potato sample was located in a microwave electromagnetic field. The diameter and height of the sample were 45 mm and 15 mm, respectively, and it consisted of solid phase, liquid water, and a gas mixture of water vapor and dry air. The microwave energy was incident on the top surface of the sample and decreased gradually with the incident depth. Fluid phases were transported under the gradients of concentration, pressure, and the additional flows caused by the deformation of the solid phase. The heat transportation included the conductive heat caused by the temperature gradient and the convective heat carried by the fluid phases. Phase change and dielectric volumetric heat were treated as the heat sources. As shown in Figure 1b, on the transport boundaries of the sample, water vapor was removed by the natural convection. Liquid water was extruded directly when the sample surface was saturated; otherwise, it was vaporized completely. Deformation was free under the stress caused by moisture loss and gas pressure. On the symmetric boundaries, no mass transport, heat transport, or deformation occurred. The following assumptions were made to simplify the problem: (1) The sample was treated as a multiphase porous medium consisting of solid phase, liquid water, and a gas mixture of water vapor and dry air. (2) The sample was treated as a compressible hyperelastic material, and the volume of the solid phase remained constant during the whole process. (3) The deformation was a quasi-state process. (4) Local thermal equilibrium existed among different phases. (5) The distributions of the temperature, moisture, and stress in the sample were initially uniform. (6) The state parameters of the environment (i.e., temperature, pressure, humidity) remained constant during the whole process. (7) All gas phases obeyed the ideal gas law and shared the same volume. (8) Phase change between the liquid water and water vapor was a non-equilibrium process. (9) Bound water and gravity were neglected.

Descriptions of the Variables
The theoretical model was developed on the basis of a representative elementary volume, which was the summation of the volume of all phases, The total porosity is the volumetric fraction occupied by the fluid phases, namely, the summation of the porosity of the liquid water and that of the gas mixture, Moreover, the total porosity at any time can be determined by the constant volume of the solid phase, The saturation of the liquid water and the gas mixture was defined as the volumetric fraction occupied by the liquid water or the gas mixture in the pore regions, Additionally, the mathematic expression of the pressure of the water vapor and dry air can be obtained on the basis of the ideal gas law, The moisture content on the dry basis can be calculated by

Conservation Equations
The total mass fluxes of the fluid phases were attributed to three parts: (1) The fluxes due to the concentration differences of different locations. (2) The fluxes flowing along the reverse gas pressure gradients. (3) The fluxes flowing with respect to the ground frame of the reference due to the solid movement. Therefore, the mass conservation equations of liquid water and water vapor were expressed as

The Non-Equilibrium Description of Phase Change
The process of the phase change represents the conversion of the liquid water and the water vapor, which causes changes in the internal gas pressure. It is described by a non-equilibrium approach [18], which obtains the evaporation rate through the difference between the calculated vapor pressure and the equilibrium vapor pressure, where the water activity and saturation pressure of the water vapor are expressed as functions of moisture and temperature [18,51],

Initial and Boundary Conditions
The initial conditions of the mass conservation equations are given in Table 1. On boundaries 1 and 2, there was no fluid flow. On boundaries 3 and 4, the removed moisture consisted of the flows caused by the concentration differences, the pressure gradients, and the solid deformation; thus, The gas pressure was obtained by solving the conservation equation of the gas mixture,

Initial and Boundary Conditions
The initial condition for the conservation equation of the gas mixture is presented in Table 1. On boundaries 1 and 2, the pressure was zero, and on boundaries 3 and 4, the pressure was equal to the ambient pressure, namely,

Conservation Equations
A single equation is enough for solving the local temperature, since thermal equilibrium is assumed to exist at one location among different phases; thus, The effective thermal properties used in the equation of energy conservation can be obtained by the weighted average of individual component,

Heat Fluxes and Sources
There were two driving forces within the sample to push the heat flow, namely, the convective heat carried by the transported fluid phases and the conductive heat due to the temperature gradients at different locations. During the thermal process, the main heat sources included the latent heat of phase change and the volumetric dielectric heat generated by the microwave interacting with the dielectric material. Additionally, the microwave power was assumed to exponentially decay along the depth direction on the basis of Lambert's law; thus, where the attenuation constant can be calculated as

Initial and Boundary Conditions
The initial condition of the energy conservation equation is given in Table 1. In boundaries 1 and 2 there was no heat flow. In boundaries 3 and 4, four forms of heat fluxes occurred, namely, the heat exchanging with the environment, the heat carried by fluid phases, the heat used to phase change, and the heat radiating to the environment,

Solid Momentum Balance
In this study, the sample deformed more than 10% of the initial volume, and thus a large deformation analysis had to be taken into consideration [19,42,43]. The driving forces of the deformation were the moisture loss and the gas pressure gradient. Therefore, the total deformation gradient tensor was characterized by a multiplicative split of the elastic deformation gradient tensor and the moisture loss deformation gradient tensor [42,43], In addition, the relation of the determinants of the total deformation gradient tensor, the elastic deformation gradient tensor, and the moisture loss deformation gradient tensor can be presented as The total internal force was zero due to the quasi-static process assumption. Therefore, the sum of the internal stress and the body force resulting from the gas pressure gradient was zero in terms of the solid momentum balance. Moreover, the Cauchy stress is related to the second Piola-Kirchoff (PK2) stress; thus, PK2 stress tensor is an energy conjugate to the Green-Lagrange strain tensor [19,42,43,45], Thus, the PK2 stress tensor and the Green-Lagrange tensor can be related via this strain energy density on the basis of the second law of thermodynamics,

Constitutive Law
In this study, the stress-strain relationship of the sample behaved as a hyperelastic material, and thus a Neo-Hooken constitutive model was chosen to obtain the strain energy density [42,43], The deformation due to the water loss is based on the assumption that the volume of the removed liquid water is occupied by the solid phase; thus,

Initial and Boundary Conditions
The initial condition of the solid momentum balance was zero. Moreover, there was no displacement on boundaries 1 and 2, while on boundaries 3 and 4, it was free to deform during the IMW thermal process.

Input Parameters
The constant input parameters were summarized and given in Table 1. The input parameters related to variables or empirical formula are discussed and presented in this subsection.

Transport Properties
The molecular diffusivity of water vapor in a porous medium is related to the gas pressure, temperature, porosity, tortuosity, and the saturation [17]. The transport of liquid water due to the capillarity can be written in terms of capillary diffusion [42,43], and the diffusion coefficient is a function of moisture content and temperature [52]; thus, When multiphase flowed through a porous medium, intrinsic permeability and relative permeability were introduced. The intrinsic permeability of liquid water in potato tissue was 1 × 10 −15 m 2 , and the intrinsic permeability of gas phase was calculated from that of liquid water via the Klinkenberg correction factor [43]; thus, The relative permeability of the liquid water and the gas phase were the functions of water saturation [43], Moreover, the permeability of the fluid phase was affected by the porosity during the thermal process. Therefore, the total permeability was corrected by a factor given by the Kozeny-Carman equation [43], In a microwave oven, there is typically a gentle natural convective air flow, as reported in the literature. Thus, we assumed a laminar flow of the air in the microwave oven, and the heat transfer coefficient was calculated as [25] h t = 0.0296 · R 0.5 e · P 0.33 The Chilton-Colburn analogy was used to calculate the convective mass transfer coefficient [43], The viscosities of the liquid water and gas mixture can be presented as functions of temperature, The thermal transport properties of liquid water used in this model were also expressed as the functions of temperature [42],

Dielectric Properties
The Landau-Lifshitz-Looyenga equation (LLLE) has been used to predict the dielectric properties of some food materials for a range of microwave frequencies with sufficient accuracy [43], where the dielectric constants and dielectric loss factors of the solid and the gas phases are constants shown in Table 1, and those of liquid water are expressed as functions of temperature,

Mechanical Properties
In order to obtain the strain energy density, one requires the Lamé's constant and shear modulus [53,54], The elastic modulus can be measured as a function of moisture content, and is given as [42] The potato sample exhibited a soft and rubber-like behavior when its temperature was above the glass transition temperature. Namely, the ideal shrinkage occurred during this period of time, and the Poisson's ratio was 0.49 (in order to have a better convergence). When the temperature was below the glass transition temperature, the sample exhibited a hard and glass-like behavior, and the Poisson's ratio was 0.3 [42,46], namely, The glass transition temperature of potato starch can be given as a function of moisture content [46],

Overview of the Validation Experiments
In order to verify the correctness of the theoretical model, we performed validation experiments using a domestic microwave oven. A fresh potato was cut into pieces with a cylindrical mold; the diameter and height of the sample are shown in Figure 2d. Before being heated in the microwave oven, the potato cylinders were boiled for five minutes in hot water for preventing enzymatic browning, and then cooled to room temperature. When the IMW process was performed, a single potato cylindrical sample was placed in the middle of a domestic flat plate microwave oven, as shown in Figure 2a,b, and then the IMW was loaded manually. The initial and final states of the potato sample are shown in Figure 2c,e. An obvious shrinkage could be observed.

Description of the Validation Experiments
The potato sample was heated by a domestic microwave oven (Midea, EG 720FC8-NR1, Guangzhou, China). The mass change of the sample was measured by an electronic balance (YINGHENG, Zhejiang, China) with an accuracy of 0.01 g. The variations of the radial and axial dimensions were measured by a vernier caliper (LINKS, HMCT GROUP, China) with an accuracy of 0.01 mm. The temperature at the surface center of the cylindrical potato was measured by a hand-held infrared thermometer (Fluke, F59C, Fortive, Anhui, China) with an accuracy of ±2 • C.
One intermittent cycle was 100 s with 20 s power-on heating stage and 80 s power-off tempering stage. After the heating stage of each cycle, the sample was moved from the microwave oven to measure and record the mass, top center temperature, and the diameter and height. The final experimental values were the mean of three replicates.

Statistical Analysis and Measurement Uncertainty
The statistical analysis and measurement uncertainty analysis were conducted by using Excel 2007. The standard deviation of the three replicate experiments and the accuracy of the measuring tools were treated as the main sources of the measurement uncertainty. Thus, the measurement uncertainty was described by Equations (51)- (54). The measurement uncertainty is presented as error bars on the curves of the experimental data in the next section.
In addition, the commonly used statistical indices (i.e., the coefficient of determination R 2 and the root mean square error RMSE) were used to evaluate the consistency of the simulation results of the theoretical models with the experimental data. Thus, the coefficient of determination and the root mean square error were described as

Solution Procedure
The multiphase transport and deformation in the porous dielectric materials during the thermal process were solved by the commercial finite element software COMSOL Multiphysics. The variables solved in different physic fields were fully coupled. MUMPS direct solver was used for the time-dependent study. A computer with a 2.10 GHz processor and 64.0 GB RAM was used to run the software.
In order to obtain the rapid and accurate simulation results, one should choose a suitable mesh size. In the finite element numerical solution of the microwave field, the maximum mesh size in the dielectric materials had a relationship with the free space wavelength and the dielectric constant [55] as Thus, we set the maximum mesh size obtained by Equation (57) as a start and refined the mesh step by step. In addition, the mesh sizes on the transport boundaries were refined due to the transport of mass and heat. As shown in Figure 3a, the average moisture content curves obtained by different mesh sizes are presented. From Figure 3a, the values of average moisture content at the same time point tended to be stable with the increase in the solving degrees of freedom. In order to save computer resources and calculation time, we selected the solving degree of freedom of 84,249 marked by the dash line in this study. The corresponding maximum mesh size was 0.1125 mm on the transport boundaries and 0.225 mm in the remaining area, shown in Figure 3b.   The average moisture content obtained by Model 1 decreased from 4.374 to 0.137, being located within the range consisting of the measuring values and the measurement uncertainty. It matched better with the experimental data with a higher coefficient of determination of R 2 = 0.987 and a lower root mean square error of RMSE = 0.180. The results of the average moisture content obtained by Model 2 decreased from 4.374 to 1.001, which was higher than the experimental data with a lower coefficient of determination of R 2 = 0.739 and a higher root mean square error of RMSE = 0.678. This was because during the thermal process, the deformation caused an additional flow of the fluid phases. Thus, the results of Model 1 were more consistent with the experimental data, namely, it was more in line with the actual situation. A similar conclusion has also been drawn by Gulati et al. [43].
In order to improve the reliability of the model developed in this study, we applied Model 1 to simulate the intermittent microwave convective drying of apple slice studied in [18]. The average moisture content obtained by Model 1 was compared with the experimental results obtained in [18], as shown in Figure 4b. A good consistence of the two results can be observed. Thus, it indicates that the model developed in our study could be extended to the application of other microwave thermal processes and could gain accurate simulation results.

Deformation Displacements
As shown in Figure 5a,b, the simulation results of the radial and axial displacements are presented and compared with the experimental data. The error bars on the curves of the experimental data represent the values of the measurement uncertainty. The negative values indicated the shrinkage of the sample. The radial and axial displacements obtained by Model 1 were −7.04 mm and −3.39 mm, respectively; the values of the coefficient of determination of the radial and axial displacements R 2 = 0.963 and R 2 = 0.926, respectively; and the values of the root mean square error of the radial and axial displacements were RMSE = 0.506 mm and RMSE = 0.335 mm, respectively. Most of the simulation results were located within the range consisting of the measuring values and the measurement uncertainty. Therefore, the simulation results of the radial and axial displacements matched well with the experimental data. Moreover, the deformation rates are also given in Figure 5, which indicates that the deformation rate reached the global maximum around 600 s due to the collapse of cell structure caused by moisture loss [56].

Temperature
In Figure 6, the temperature results at the top surface center of the sample obtained by Model 1 and Model 2 are compared with the experimental data, and the error bars on the curves of the experimental data represent the values of the measurement uncertainty. The temperature obtained by the experiment ranged from 29.9 • C to 73.65 • C, and the values of the temperature obtained by Model 1 and Model 2 ranged from 30.0 • C to 86.17 • C and ranged from 30.0 • C to 85.34 • C, respectively. Most of the simulation results were located out of the range consisting of the measuring values and the measurement uncertainty. Moreover, the coefficients of determination and the root mean square errors of the two models were R 2 = 0.437 and R 2 = 0.345, and RMSE = 13.5 • C and RMSE = 15.4 • C, respectively. This indicates that the simulated temperature values were higher than the experimental data. However, there was temperature deviation between the simulation results and the experimental data, and the temperature trends of both models were similar to the experimental data. Thus, the models can be used to describe the temperature tendency. There were two reasons for this large temperature deviation between the simulation and experimental results. One was the inevitable measuring error of the hand-held infrared thermometer; the other reason was the usage of Lambert's law to calculate the electric filed. Lambert's law is based on the semi-infinite critical sample length assumption, and it does not consider the standing wave effects [57]. Moreover, the power absorption on the surface is always the maximum, regardless of the moisture content [18,58]. Similar overestimations of the temperature obtained by the models using Lambert's law have also been observed in previous studies, and the Maxwell's equations are more ideal in modelling the microwave heat generation and the standing wave effects [18,57]. Therefore, we used the Maxwell's equations to simulate the distribution of the electromagnetic field within the sample in our next stage of study.

The Mass Transport Mechanism
The curves in Figure 7a,b present the time evolution of the average values of the mass fluxes of different fluid phases and those under different driving forces, respectively. As shown in Figure 7a, initially, the average flux of liquid water was higher than that of water vapor by three orders of magnitude on average; then, the average flux of liquid water turned to being lower than that of water vapor by one order of magnitude on average. Thus, combining the time evolution curves of the moisture content, as shown in Figure 4a, the removing of liquid water was first the main part of the entire moisture removing, which made the moisture content decrease from 4.37 to 1.11; then, the removing of water vapor turned to being predominant, which made the moisture content decrease from 1.11 to 0.14. This is consistent with the statements in [16,17]. Therefore, during the IMW thermal process, when the moisture content was high, the moisture loss was caused by the removal of the liquid water first; then, it turned to water vapor when the moisture content was low. As shown in Figure 7b, the values of the diffusion flux driven by concentration gradients and those of the convective flux driven by pressure gradients were on the same order of magnitude, and the values of the former were slightly higher. Thus, the diffusion flux contributed more to the moisture loss. In addition, the Peclet number (Pe) defined as the ratio of the convective (mass/heat) transport to the diffusion (mass/heat) transport [59,60] was used to reflect the dominant mechanism of the (mass/heat) transport. In Figure 8, the time evolution curves of the maximum and average mass Pe and the distribution diagrams of the mass Pe at some specific time points are presented. In each cycle, the mass Pe increased to a peak value after the heating stage, then decreased to a valley value after the tempering stage. The peak values of the maximum mass Pe ranged from 0.04 to 0.20, and those of the average mass Pe ranged from 0.01 to 0.05. All the peak values were less than 1.00. Namely, the mass transport by diffusion driven by the concentration gradient was larger than that by convection driven by the gas pressure gradient. Thus, the diffusion was the dominant mechanism of the mass transport during the whole IMW thermal process, which is consistent with the conclusion drawn from Figure 7b. From the mass Pe distribution diagrams shown in Figure 8, the mass Pe was relatively large near the surface due to the higher velocity of the liquid phase.

The Heat Transport Mechanism
The curves in Figure 9a,b demonstrate the average heat fluxes of different phases and different heat transfer methods, respectively. In Figure 9a, the heat fluxes of liquid water and solid phase were on the same order of magnitude, and both were higher than those of water vapor and dry air by one order of magnitude. Thus, the liquid water and the solid phase were the principal media carrying heat. From Figure 9b, we see that both convective and diffusion heat fluxes were on the same order of magnitude, and the two values were comparable. Except for the period between 300 s and 700 s, the peak values of the average diffusion heat flux were higher than those of the average convection heat flux. The valley values of the average diffusion heat flux were always higher than those of the average convection heat flux during the whole process. In Figure 10, the time evolution curves of the maximum and average heat Pe and the distribution diagrams of the heat Pe at some specific time points are presented. In each cycle, the heat Pe increased to a peak value after the heating stage, then decreased to a valley value after the tempering stage. The peak values of the maximum heat Pe ranged from 0.14 to 23.15, and those of the average heat Pe ranged from 0.03 to 1.73. Most of the peak values of the maximum and average heat Pe were larger than 1.00, and all of the valley values were lower than 1.00. Namely, the convective heat flux was the dominant mechanism in most of the heating stage, and then the diffusion heat flux turned to being dominant in the tempering stage of all cycles. Furthermore, as shown in Figures 8 and 10, the peak values of both mass and heat Pe reached local maxima in the middle and final stages of the whole processes, respectively. In the middle stage, the increase in the peak values of Pe resulted from the increase in the velocities of liquid water and solid deformation. On the other hand, in the final stage, the increase in the peak values of Pe was caused by the increase in the velocity of water vapor.

The von Mises Stress
The von Mises stress, based on the fourth strength theory, is used to predict the quality attribute of the sample (e.g., the failure or crack of the material). The time evolution curves of the maximum and average von Mises stress are shown in Figure 11. Moreover, the von Mises stress distribution diagrams at some specific time points are presented in Figure 11. In each cycle, the von Mises stress increased to a peak value after the heating stage, then decreased to a valley value after the tempering stage. The peak values of the maximum von Mises stress ranged from 0.01 MPa to 0.36 MPa, and those of the average von Mises stress ranged from 0.00 MPa to 0.15 MPa. It revealed that the curves of the maximum and average von Mises stresses shared a similar trend. Namely, the peak values of both stresses reached the local maxima of 0.18 MPa and 0.05 MPa at 520 s, respectively; then, the peak values of both stresses reached another local maxima of 0.36 MPa and 0.15 MPa at 1420 s, respectively. In addition, the local maximum stress at 520 s appeared on the sample surface, while that at 1420 s was in the material. The reason for the first local maximum stress was the high gradients of temperature, gas pressure, water vapor concentration, and liquid water concentration [43]. On the other hand, the second local maximum stress was due to the glass transition of the material, which made the elastic modulus larger and the Poisson's ratio lower. Moreover, material failure or crack would appear at a location if the stress there was above the fracture stress. Thus, the crack or fracture most likely occurred at the locations with the local maximum stress. Therefore, the fracture most likely appeared on the sample surface at 520 s or inside the materials at 1420 s. On the other hand, the fracture stress of potato was on the order of 1 × 10 6 Pa magnitude [43,52], which was higher than the maxima von Mises stress in the sample shown in Figure 11. Therefore, in the sample of our study, crack or fracture of the product can be avoided, and drying product can be obtained in a good quality. This was validated by the final state of the potato sample shown in Figure 2e.

The Effects of the Cycle Length on the Outcomes
The power ratio or pulse ratio (PR) is a critical parameter used in the intermittent process. It is defined as the ratio of the total time length of a cycle to the heating time in the cycle or its inverse [61][62][63]. From previous studies, it was revealed that prolonging the tempering time is of benefit to the quality improvement of the products [61][62][63]. However, there is no study on the effects of the cycle length on the results of the process with the same PR. Thus, in this section, we focus on the effects of the cycle length on the results of the porous dielectric material in the IMW thermal process. The settings of different cycle lengths are shown in Table 2.  , s  2  4  8  12  16  20  24  Tempering, s  8  16  32  48  64  80  96  PR  5  5  5  5  5  5  5  Cycles  120  60  30  20 15 12 10 In Figure 12a,b, the solid lines are the time evolution curves of the average moisture content on dry basis and the maximum displacement, respectively; the dot-dashed lines are the values of the average moisture content and the maximum displacement of different settings at 1200 s. From Figure 12a,b, it can be seen that the fluctuations of the time evolution curves increased with the prolonging of the cycle length. The changing rate of the longer cycle length was also higher. The average moisture of Set 1 was the lowest, changing from 4.37 to 0.29, and that of Set 3 was the highest, changing from 4.37 to 0.51. Similarly, the maximum displacement of Set 1 was the highest, with a changing range from 0.00 mm to 7.55 mm, and that of Set 3 was the lowest, changing from 0.00 mm to 7.15 mm. At 1200 s, the final values of the average moisture content and the maximum displacement of different settings presented fluctuations with the cycle length. Namely, the final average moisture content of Set 1 was the lowest at 0.29, and that of Set 3 was the highest at 0.51. The final maximum displacement of Set 1 was the highest at 7.55 mm, and that of Set 3 was the lowest at 7.15 mm. In addition, we used the standard deviation (SD) to evaluate the gap among different simulation results, as shown in Figure 12c,d. It indicated that the gap among different settings reached the global maximum between 500 s and 600 s, which resulted from the high water loss and deformation rates at this stage.
In Figure 13a,b, the time evolution curves of the peak and valley values of the average temperature are presented, and those of the maximum von Mises stress are presented in Figure 13c,d. From Figure 13a, the peak values of the average temperature of different settings shared a similar trend, which increased firstly and then kept stable. As a whole, the peak values of the average temperature of different settings increased with the cycle length. Similarly, the valley values of the average temperature of different settings decreased with the cycle length, as shown in Figure 13b. In Figure 13c,d, the peak and valley values of the maximum von Mises stress also increased and decreased with the cycle length as a whole, respectively. In addition, both the peak and valley values of the von Mises stress both reached local maximum values around 500 s, which resulted from the high gradients of the moisture content, temperature, deformation, and pressure at this stage [43]. Therefore, from the above comparisons of the results of different cycle lengths, we found that the fluctuations of the outcomes became larger with the increase of the cycle length.

The Effects of the Frequency on the Outcomes
The differences between the radio frequency (RF) and the microwave thermal processes were all caused by their frequency difference [8]. The model developed in our study can be extended to simulate the IRF thermal process. We compare the results of the IMW and the IRF thermal processes in this section. The time evolution curves of the average moisture content, the maximum displacement, the average temperature, and the maximum von Mises stress are presented in Figure 14a-d, respectively. Moreover, the absolute differences between the results of the IMW and IRF (i.e., the values obtained by the IRF minus those obtained by the IMW) are also presented in these figures. As shown in Figure 14a, the average moisture content of the IRF was lower than that of the IMW. Namely, the average moisture content of the IRF decreased from 4.37 to 0.23, and that of the IMW decreased from 4.37 to 0.45. In Figure 14b, the maximum displacement of the IRF was higher than that of the IMW. Namely, the maximum displacement of the IRF increased from 0.00 mm to 7.66 mm, and that of the IMW increased from 0.00 mm to 7.27 mm. Therefore, if the dielectric materials need to be dried with a higher efficiency, the IRF is a better choice, while if the purpose of the thermal process is heating with less water loss and deformation, then the IMW is more ideal instead. In addition, both of the absolute differences of the average moisture content and the maximum displacement of the two processes had global extremes around 500 s, being −0.91 and 1.63 mm, respectively. This was due to the violent transport of mass and heat at this period of time.
As shown in Figure 14c,d, in the heating stage of each cycle, the average temperature and the maximum stress of the material treated by the IRF increased faster than those treated by the IMW. In the tempering stage of each cycle, the average temperature and the maximum stress of the material treated by the IRF decreased faster than those treated by the IMW. In Figure 14c, the average temperature peak values of the IRF were higher than those of the IMW. While the average temperature valley values of the IRF were higher than those of the IMW before 600 s, they then turned to be lower thereafter. The absolute differences of the average temperature between the two processes reached a global maximum of 27.22 • C at 520 s. Therefore, the temperature of the material treated by the IRF increased faster than that treated by the IMW; however, the average temperature peak values of the materials treated by both processes tended to reach a similar equilibrium and stayed constant thereafter. In Figure 14d, the peak values of the maximum stress of the IRF was higher than those of the IMW, except for the period between 420 s and 700 s, while the valley values of the IRF were higher before 400 s, and then the values were lower than those of the IMW thereafter. The maximum stress peak values of both processes reached local maxima in the middle stage and then decreased to the local minima, and the local maximum of the IRF reached earlier than that of the IMW. Namely, the stress peak values of the material treated by the IRF reached a local maximum of 0.19 MPa at 420 s, then decreased to a local minimum of 0.12 MPa at 620 s, and those of the material treated by the IMW reached a local maximum of 0.18 MPa at 520 s, then decreased to a local minimum of 0.12 MPa at 720 s. The stress valley values of both processes reached a global maximum of 0.08 MPa at 400 s. Therefore, the maximum stress in the material treated by the IRF was similar to that treated by the IMW in the aspect of amplitudes and trends, both were lower than the fracture stress. The difference between the two processes was that the local maximum of the IRF arrived earlier than that of the IMW, which was due to the deeper electromagnetic energy penetration of the IRF.

Conclusions
A 2D axi-symmetric physically based theoretical model was developed, and the electromagnetic energy, multiphase transport, phase change, large deformation, and glass transition were taken into consideration. The model was validated by an IMW experiment, as well as the experimental results of other researchers. Therefore, the model can be used in various processes (e.g., microwave thermal processes, radio frequency thermal processes, and even other thermal processes) of different porous media. When the model is employed in a specific industrial application, only the physical parameters of the corresponding materials and the environment need to be set. It can uncover the mechanisms of the mass and heat transport intuitively and can be used to optimize the process parameters in an efficient and economical way on the basis of the simulation results of temperature, moisture content, and deformation.
From our simulation results and discussions, it can be concluded that the mass was mainly carried by liquid water when the moisture content was high, and turned to being carried by water vapor when the moisture content became low. The heat was mainly carried by liquid water and solid. By comparing the average fluxes and Pe numbers of the mass and heat, we found that the diffusion was the dominant mechanism of the mass transport during the whole process. For the heat transport, the diffusion was the dominant mechanism during the tempering stage, while the convection dominated the heat transport near the surface areas of the material during the heating stage. In our study, throughout the whole process, the von Mises stress was lower than the fracture stress. The von Mises stress reached a local maximum on the top surface of material in the middle stage due to the high gradients of the moisture content, deformation, and temperature, and it reached another local maximum inside the material at the final stage due to the glass transition of the material.
For the intermittent thermal process, different cycle lengths with the same PR affected the outcomes. In our study, the setting of 10 s/cycle led to the outcomes with the lowest moisture content and the largest deformation, the setting of 30 s/cycle led to the outcomes with the highest moisture content and the least deformation, and the outcomes of other settings lay between them. Moreover, the fluctuations of the outcomes of the moisture content, deformation, temperature, and von Mises stress became larger with the increasing of the cycle length. Thus, the material treated by a longer cycle length tended to trigger the phenomena of overheat and fracture.
The model in our study can be used to simulate both the IMW and IRF thermal processes. By comparing the results of the two thermal processes, we found that the IRF process tended to obtain the lower moisture content and the larger deformation. Moreover, the outcomes of the temperature and von Mises stress in the two processes were similar. The only difference was that the IRF reached the local maxima earlier than the IMW due to the deeper penetration of the electromagnetic energy. Thus, for the drying process, the IRF was found to be more suitable due to its lower moisture content outcome, while for the heating process requiring less deformation, the IMW can be said to be more ideal.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to unpublished results.

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.