Modeling of Dielectric Heating within Lyophilization Process

A process of lyophilization of paper books is modeled. The process of drying is controlled by a dielectric heating system. From the physical viewpoint, the task represents a 2D coupled problem described by two partial differential equations for the electric and temperature fields. The material parameters are supposed to be temperature-dependent functions. The continuous mathematical model is solved numerically. The methodology is illustrated with some examples whose results are discussed.


Introduction
From the technological viewpoint, sublimation is cleaning or drying method which represents direct phase transition from the solid to the gas-phase without melting, passing through the liquid phase and boiling.Nearly every given material can sublime at certain conditions.These conditions are usually set from a proper phase diagram.There are two different types of sublimation-at the normal pressure and vacuum sublimation.
Sublimation at the normal (or atmospheric) pressure is usual for volatile matters, such as iodine or solid carbon dioxide (dry ice).Sometimes it happens that the room temperature is not high enough to cause this effect even by very volatile matters.In that case careful warming is needed.This type of sublimation acts for cleaning technologies and processes only (associated with the backward condensation).
Vacuum sublimation is faster, gentler, often thriftier, and technologically more controllable.It can act for drying processes too.This is the case of the discussed technology called lyophilization.Lyophilization is a dehydration process at the conditions of low temperature and reduced pressure.In fact, the process works by freezing the proper substance and is also known as freeze-drying.The mentioned conditions enable the substance to sublime.The process is heavy on energy and demands special apparatus called freeze-drier [1,2].This paper deals with the possibility of potentially new technology based on the freeze-driers supplemented with a dielectric heating system for controlling the process inside.The inner space of the freeze-drier (Figure 1) is connected with a vacuum pump and an ice condenser, where the newly emerged vapor is concentrated.The process of lyophilization is commonly divided into the following stages.
The first one is freezing at normal pressure and it is done using a freezing machine.The final freezing temperature is ordinarily set somewhere in the range from −15 to −50 ∘ C because it is necessary to cool the material below its triple point.The rate of the freezing is very important, too.Larger crystals are produced by slower freezing and are easier to freeze-dry.However, larger crystals have more destructive effect on the material than smaller ones.For example, in the case of drying old paper books, this is the decisive reason to do the freezing rapidly, in order to lower the temperature of material below its eutectic point quickly.Therefore, the freezing stage is the most critical in the whole process.
During the next stage-primary drying-the pressure is lowered to the range between 100 and 500 Pa.This evacuation is held at the proper temperature for the specific material.For the beginning of the sublimation of ice (water) a sufficient amount of heat has to be supplied to the material.In this first drying stage, about 95% of the water content in the material is sublimated.This stage is usually slow because not too much heat is allowed to be added.The pressure inside is controlled by the vacuum pump.The reduced pressure speeds up the process.Furthermore, the ice condenser gives a surface for the water vapor to resolidify on.It does not have any other meaning in the process, but it prevents water vapor from reaching the vacuum pump.The temperature in the condenser is not necessarily lower than the temperature in the freeze-drier.This stage is the most important for the examples of the paper and will be discussed later.Therefore, it should be noted that the heat is brought mainly by conduction and radiation.Due to the reduced pressure, the effect of convection is negligible [3,4].
Finally, the task of the secondary drying stage is to remove liquid water particles-the rest after the primary drying.This stage is driven by sorption isotherms of the given material.The temperature is raised and can be even above 0 ∘ C.However, the pressure in the drier remains usually constant, so the process of the vacuum sublimation is transformed just a little.Sometimes the pressure is even lowered in the range of fractions of a Pascal.But it is not a rule.In some cases there is a benefit from increased pressure as well.When the process is complete, the inner pressure is broken with an inert gas (e.g., argon, nitrogen) or filtered dry air.After this last process, the final residual water content in the material is around 1% to 4% [1][2][3][4][5][6].
Lyophilization is a very often used process nowadays but usually connected with resistance heating.The idea of utilization of the dielectric heating is relatively novel.The reason for the utilization is the different mechanism for drying with dielectric energy.With internal heat generation, most of the moisture is vaporized before leaving the material.The drying is very rapid even without the need to overheat the space.In general, the dielectric heating can be found in numerous systems in the lumber, furniture, textile, paper, food, tire, and ceramic industries [7][8][9][10] nowadays.Focused on the drying of the books, the paper industry is the most interesting one.The dielectric heating is used to dry printing inks, adhesives, and coating materials on paper and to dry the paper itself.This is possible to combine with hot air, infrared, or other heating media to achieve optimum results.For papermaking, machines are in use with 50 to 500 kW [11,12] of radio frequency at 13.56 or 27.12 MHz [13].The radio frequency is able to overcome the most common problems encountered at the dry end of the process-low efficiency due to moisture distribution through the paper thickness, uneven moisture distribution across the width of the web, temporary unevenness, or streaks of moisture due to some failure of the equipment.
The advantages of this technology could be found in the improved general controlling of the process.Potential energy savings achievable from such a system due to the speed of drying are of great interest as well [14][15][16].The technology seems to be more efficient and, if well controlled, nondestructive.The drying can be generally done at low ambient temperatures.This paper aims to verify these theoretical advantages.Since the primary drying is the heaviest stage on energy, the paper is focused on that, as already mentioned above.All procedures and theoretical issues are quoted and discussed in the very next part.

Mathematical Models
In fact, the stated coupled problem of dielectric heating consists of two physical fields-on the one hand of electric field, on the other hand of temperature field.For more transparent conception of the theoretical issues and procedures of the calculations in future examples both fields are discussed separately.

Electric Field.
For the description of electric field at the dielectric heating, it is considered that the heated material is situated between two electrodes of a parallel-plate capacitor.The geometry of the capacitor is characterized by the surface of the electrodes  and their constant distance .If the deformation of the field at the margins of both electrodes is neglected, a homogenous electric field can be considered.In the next parts few examples about the freeze-drying of books will be shown.The book paper is a typical dielectric material with its temperature-dependent parameters, such as tan() (discussed in the very next paragraph), relative permittivity   , thermal conductivity , and the product of material density and its specific heat capacity at constant pressure  ⋅   .
Every real capacitor is characterized by its losses and can be considered as a parallel connection of an ideal resistor (with its resistance ) and an ideal capacitor (with capacitance ).The real capacitor is connected to the source of harmonic alternating voltage.Electric current flowing through the equivalent circuit is divided into the resistance component   and the component of the ideal capacitance   .The angle  which denotes the relation of the voltage  to the general electric current , here both in the form of phasors, is determined as where  stands for the angular velocity which depends on the nominal frequency of the dielectric heating : in (1) can be easily substituted for well-known relation describing the capacitance of the parallel-plate capacitor which leads to the formula where is the vacuum permittivity.Anyway, the resistance of the dielectric material  can be also determined as where  denotes the electrical conductivity of the material.Finally, the last substitution of  in (3) for ( 5) transforms (3) into which gives The volumetric thermal losses   can be obtained from the differential Ohm's law as where  stands for the electric field strength.The electrical conductivity  can be easily substituted for (7).It leads to another simple expression of the volumetric losses: Symbol   is just another signage of the same meaning.It is mentioned due to the conventions in heat transfer equation, as discussed later.
The mentioned expressions and equations describe some relations of the given material parameters to the functions of the electric field, such as the electric displacement  or the electric field strength .However, the distribution of electric field in space is correctly described assuming where  is the current density.Regarding the differential Ohm's law again, (10) is transcribed to the form which leads to the following partial differential equation: where (, ) denotes the electric potential as a function of position between the electrodes.It is the final important equation for the determination and recalculation of the electric field after a fixed time step  step .This step will be defined for the variable material parameters affected by the changing temperature field.The nominal voltage   of the dielectric heating is expressed as a difference between the potentials at both electrodes: The nominal voltage   will be given.It means that the potentials define the boundary conditions.Since the potential (0, ) can be set as the ground potential, the significance of (, ) is clear from (13).
2.2.Temperature Field.The most general description of the temperature field is Fourier-Kirchhoff heat transfer equation which is considered in the form where  stands for the density of the dielectric material,   is the specific heat capacity at constant pressure,  denotes the temperature,  is the time, v stands for the vector of velocity, and  denotes the thermal conductivity.  is already known from (9) [17,18].Since any part of the system does not execute any movement, the velocity v = 0 and ( 14) can be simplified and mentioned as which is still a parabolic partial differential equation.It is needed to be supplemented with initial and boundary conditions.The initial condition sets the initial temperature (, , ) = (, , 0) in the freeze-drier at the starting point of the process of the dielectric heating as The boundary conditions of the temperature (, , ) are more complicated.They have to be given in the form of the following equations to set the temperatures (0, , ) and (, , ) at both electrodes as functions of time: where  left and  right stand for the heat transfer coefficients of the first and the second electrode,  outleft and  outright denote the values of temperature in the near surround of both electrodes,  is the emissivity of the electrodes, and is the so-called Stefan-Boltzmann constant.

Numerical Solution
The principal computations are carried out considering the task as a 2D coupled problem.Regarding the very next part of this paper, the distribution of the temperature in the Cartesian -direction of the book with constant  coordinate somewhere in the central area (half of the height approximately) could be described by a 1D simulation because the distribution is not affected by shapes of the top and the bottom.This property was utilized for checking of the correctness and results of the 2D model.These checking 1D computations were performed using the code Wolfram Mathematica, where the electric field was solved by the shooting method and the temperature field by the numerical method of lines, in comparison with the 2D simulation performed using the code Agros2D, where the problem was solved by the finite element method.Whereas the finite element method (FEM) is well-known and its discussion is not the subject of this paper the shooting method and the numerical method of lines may not be well-known and it is the reason for a short introduction of them in the next paragraphs.
The shooting method works by considering the boundary conditions as a multivariate function of initial conditions at some point, reducing the boundary value problem to finding the initial conditions that give a root.The advantage of the shooting method is that it takes advantage of the speed and adaptivity of methods for initial value problems.The disadvantage of the method is that it is not as robust as finite difference or collocation methods: some initial value problems with growing modes are inherently unstable.The shooting method looks for constant initial conditions.
The numerical method of lines is a technique for solving partial differential equations by discretizing them in many dimensions, and then integrating the semidiscrete problem as a system of ordinary differential or differential algebraic equations.A significant advantage of the method is that it allows the solution taking advantage of the sophisticated general-purpose methods and software that is developed for numerically integrating given differential equations.For the partial differential equations to which the method of lines is applicable, the method typically proves to be quite efficient.It is necessary that the partial differential problem be well posed as an initial value (or Cauchy) problem in at least one dimension, since the used integrators of the ordinary differential equations are initial problem solvers.This rules out purely elliptic equations but leaves a large class of evolution equations that can be solved quite efficiently, such as parabolic equations (e.g., temperature field).The problem is discretized with respect to variables using second-order finite differences.Even though finite difference discretizations are the most common, there is certainly no requirement that discretizations for the method of lines be done with finite differences; finite volume or even finite element discretizations can be used.The Dirichlet boundary conditions can easily be handled by a simple defining of the proper value.An alternative option is to differentiate the boundary condition with respect to time.The Neumann boundary condition is a little more difficult.It can be handled by an equivalent symmetry to itself with second-order differences [19,20].
It has been already mentioned that all of the calculations were performed within the fixed time steps.The 2D temperature field was computed within these time steps, whereas the electric field was recalculated in between.It means that it was set a new distribution of the electric field at the beginning of each step (based on the growing temperature and the changing temperature-dependent material parameters).These steps represent the coupling strategy of the model.Their chosen significance is stated among input parameters of illustrative examples.These parameters and the coupling strategy were implemented by writing a Python script within the code Agros2D.

Illustrative Examples
Three different examples are presented to investigate the potential advantages of the introduced technology.In the given arrangement the presented situations differ in the viewpoint of their symmetry from each other.The implemented asymmetries are simulated at the boundaries (electrodes, in fact) as variations of the heat transfer coefficients  left and  right considering constant temperatures  outleft ,  outright and the emissivity .
The essential input parameters are set as follows: the distance of the electrodes  = 0.05 m, the height of the book ℎ = 0.40 m, the nominal voltage   = 800 V, the working frequency  = 27 MHz, the temperature  outleft =  outright = −7 ∘ C, the emissivity  = 0.8, and the initial temperature  0 = −20 ∘ C; the chosen computational fixed time step  step = 10 s stands for a compromise between the accuracy and time demands of the computation.The accuracy of this 2D model was checked using the code Wolfram Mathematica, as already mentioned.The following material parameters, such as tan(), the relative permittivity of paper   , the thermal conductivity , and the product of the paper density and its specific heat capacity at constant pressure  ⋅   , are functions of temperature.Their temperature-dependent characteristics are shown in Figures 2, 3, 4, and 5.
All characteristics (Figures 2-5) depict estimations of the material parameters of book paper, but they are based on real values of the given material measured at normal conditions.
The point of the phase transition (sublimation) is represented by the temperature  subl = −7 ∘ C. It can be easily seen in the shown diagrams (Figures 2-5).Figure 5 illustrates the modeled behavior of the product  ⋅   at the point  subl , which denotes the amount of volumetric energy needed for the phase transition.It is possible to set this exact model considering that the emerged vapor is sufficiently and rapidly  concentrated in the ice condenser without resolidifying in the freeze-drier.
The following examples stand for the situation of symmetry, weak asymmetry, and strong asymmetry.

Example 1: Symmetrical Case.
In the case of initial symmetry both of the heat transfer coefficients are equal and set as     The mentioned examples and their consequences result from the real input parameters even in the case of the heat transfer coefficients.The asymmetries set in the second and third examples do not have such an impact as expected.The impact can be observed at the beginning of the drying process (comparing Figures 6, 10, and 14).Later, the impact of the boundary asymmetries is neglected.The rate of the neglecting effect is very transparently noticed in Figure 18.The effect disappears after 200 s, approximately.In Figures 6,  10, and 14 the effect of the dielectric heating is observable.Since the central area of the book is strongly frozen, the temperature starts growing at the electrodes.In Figures 7,8,9,11,12,13,15,16,and 17, there is a very considerable trend of flattening of the ()-characteristics by neglecting the discussed asymmetrical boundary conditions.

Conclusion
The paper presents the complete mathematical model of the dielectric heating within lyophilization process.Using the proper solver and its numerical analysis, the 2D model does not have to be solved under any simplifying assumptions.However, a variety of shapes of the book top and bottom boundaries were investigated without detecting any significantly different temperature distribution close to them (compared to the distribution in the central area).From this viewpoint, the 2D simulation is just heavy on time and does not take any additional advantage.It means that 1D model should be preferred in cases of simple book geometry.The dielectric heating turns out quite robust as a system for the lyophilization processes.It means two consequences.Firstly, initial and boundary asymmetries do not have any significant negative effect on typical smoothness of this process.Secondly, the technology can be heavy on energy.The most energy demanding moment in the process is overcoming the phase transient point, in fact the sublimation temperature (noticeable from the proper figures).But it is not the most important disadvantage.Using it for special purposes or laboratory applications, the technology can be still refundable.Nevertheless, costs of a special power supplying source are the most significant disadvantage preventing from a wider or standard industrial use.The resulting parameters

Figures 6 , 7 , 8 ,
Figures 6,7,8,and 9, depict the distribution of temperature along the -direction between the electrodes at four points in time  after the beginning of the heating process.The Cartesian  coordinate remains constant and is chosen as  = 0.20 m (half of the book height), but the results at the  coordinate close to the top or the bottom (even in case of different shapes) are equal or almost identical.