Modelling and Simulation of an Ohmic Heating Process

Ohmic Heating (OH) is one of the emerging thermal technologies used in food processing which can produce rapid and uniform heating with close to 100% energy transfer efficiency. Although mathematical modelling for OH processes has been studied by many researchers in recent years, systematic simulations of OH have not been developed for model-based control of the processes. In this paper, mathematical model for a Colinear Ohmic Heater is presented, analyzed, and studied based on the selected configuration. A numerical solution for the mathematical equations has been defined and pro-posed. MATLAB/Simulink model is hence developed and validated against the available data. Simulation results have shown that MATLAB/Simulink model can produce robust outputs at low computational costs with an accuracy of up to 99.6% in comparison to the analytical solution. This model can be used in further studies for analysis of the OH processes and development of advanced controllers.


Introduction
Heating is an essential step in food processing, it is usually used in sterilization, cooking, and pre-treatment processes of different types of foods: solid, liquid, and solid-liquid mixtures. Generally, microbiological safety and enhanced shelf life of food are the main goals of the heating process [1] [2].
Conduction, convection, and radiation are the conventional heating methods (CH) that have been commonly used in the food industry, which are based on heat transfer mechanisms, where thermal energy is transferred from one object to the other [3]. As a result of their dependency on the heat transfer mechanisms, these conventional heating methods are limited by the heat transfer rate from an external medium to the food and by the thermal conductivity of the food, which might result in over-processed products due to the lengthy processing time required to reach the target temperature, unwanted temperature peaks, and destroyed product quality [3] [4].
In order to avoid the aforementioned drawbacks and to meet the increased demand of consumers for safe, healthy, and high-quality foods items, several thermal and non-thermal technologies were investigated with the objective to meet these demands [5] [6]. Goullieux and Pain [1] pointed out that food fluids are considerably enhanced when volumetric heating methods such as Ohmic Heating (OH) are used in their heat treatment. Sastry S. [7] classified the Ohmic Heating as a Moderate Electric Field (MEF) process in which the applied electric field is less than or equal to 1000 V/cm which is considerably lower than the field strength used in the high voltage Pulsed Electric Fields (PEF) technology.
The research on the use of Ohmic Heating in food processing started by the end of the 19 th century when electricity became commercially available. Ohmic Heating (OH) or Joule heating is based on the idea of heating food through dissipating electrical energy into heat by passing an alternating current through them [8]. The first industrial use of OH was in 1920 when it was used for pasteurizing milk using carbon electrodes (220 V) and an alternating current (60 Hz).
Following that the process got approved in six US states [9]. Sarang and others [10] stated that uniform and rapid heating which results in less thermal damage in comparison with other conventional heating methods can be achieved using Ohmic Heating. Cho W. I. and others [11] showed that OH can achieve the target heating temperature of a fermented red pepper paste in 46.21 seconds in comparison with the conventional methods which took 91 minutes. Moreover, OH provided more uniform heating and produced a better microbial inactivation effect than the CH methods. In comparison with other electro-heating methods (Microwave and Radio frequency) the rapid and uniform heating resulted from the Ohmic Heating can be achieved with a lower capital cost [12], as the Ohmic Heating process has close to 100% energy transfer efficiency, while Microwave heating has 65% efficiency at its best [13].
Kamonpatana and others [4] reported that Ohmic Heating is acceptable for commercial production and that it can help to produce high-quality products. In addition to sterilization, pasteurization, and pre-treatment Ohmic Heating has a large number of potential applications such as blanching, evaporation, dehydration, fermentation, extraction, cooking, thawing, and diffusion [4] [10].
Ohmic Heating is one of the emerging thermal technologies that could be better developed, understood and validated through its mathematical model [14]. A mathematical model of the Ohmic Heating thermal process would help to ensure a completely safe ohmically cooked products, quantify heat losses, and identify the influence of the process variables [12]. A significant number of researchers introduced various mathematical models and simulations of the OH process for different types of food products (solid, liquid, and solid-liquid mixtures) and compared their modelling with the process experimental results, for instance, Shim J.S. and others [13] modelled the thermal behavior of a multiphase food products with different electrical conductivities under Ohmic Heating using computational fluid dynamics (CFD) codes, Guo W. and others [15] modelled a 3-dimensional OH for high frequencies based on Maxwell's equations of the electric field, while De Alwis and Fryer [16] used Laplace's equation which includes Joule's law and Ohm's law to model the process.
Despite all the efforts in developing a mathematical model for OH processes, very rarely models are developed and used in systematic simulations suitable for model-based control of the processes. This work aims to fill the gap with MATLAB/Simulink model being developed for better analysis, optimization, and advanced control of typical OH processes.
In this paper, the Ohmic Heating process mathematical equations are presented, analyzed, and studied based on the selected Colinear Ohmic Heater configuration. A numerical solution for solving the partial differential equations is defined and proposed. MATLAB/Simulink models are hence developed and validated against the available data in the literature.

Fundamental Ohmic Heating Mathematical Equations
First, the fundamental concept of the ohmic heating process depends on Ohm's law which states that when a current flow through a conductor (e.g. food product), the internally generated heat within the conductor is proportional to the resistance (1/Electrical Conductivity) of the conductor and the square of the applied electric intensity (Electrical Field Strength). Therefore, from the above definition, the OH process can be decomposed into two main sub-processes, the Electric process (electric field) which is the result of the applied voltage across the heater electrodes, and the process of the Internally generated thermal energy (thermal field) which is the result of the flowing current through the food product. One can observe that the two processes are strongly linked to each other, for instance, to increase the internally generated heat the strength of the applied electric field should be increased, and as the food product temperature raises its electrical resistance decreases (electrical conductivity increases), hence, more heat is generated within it with the same applied electric field strength value.
Generally, the electric field distribution within the Ohmic Heater cell can be calculated by solving Laplace's equation (Equation (1)), which states that the acceleration (the second gradient) of the electric field multiplied by the product electrical conductivity within the OH cell is equal to zero [16]. where V is the applied voltage across the electrodes and σ(T) is the food electrical conductivity as a function of its temperature. Once the electrical field is known the temperature distribution inside the food product can be calculated solving the following heat equations (Equations (2) and (3)). Equation (2) describes the relationship between the food thermo-physical properties, the internally generated energy source, and the distribution of the generated temperature, while Equation (3) defines the relationship between the internal energy source, the applied electric field, and the product electrical conductivity [17].
where ρ is the food density (kg/m 3 ), C p is its specific heat (J/kg ˚C), T is the temperature (˚C), t is the time, K is its thermal conductivity (W/m•˚C), and Q is the internal energy source (W/m 3 ) that is generated by the applied electric field (E) as illustrated in Equation (3).
As the food product temperature raises due to the applied electric field and the internally generated heat source, its electrical conductivity increases linearly as a function of its temperature, Equation (4) describes that relationship.
where σ o is the electrical conductivity at the reference temperature (initial temperature) T o , k o is the temperature coefficient (1⁄˚C), and T is the current temperature, as seen above the temperature coefficient defines how much the electrical conductivity would increase or decrease with the temperature change.
The aforementioned equations (Equations (1) and (2)) are of the partial differential type. In order to solve the OH partial differential equations which the numerical solution is obtained such as finite difference method, finite element method, and finite volume method, the selection of the numerical solution method depends mainly on the type of the equations to be solved and the system to be studied. As demonstrated in the above equations (Equations (1) to (4)) and in order to analyze the thermal behavior of the OH process simultaneous solution of the coupled electric and thermal fields should be calculated. At the first step of the solution and to solve the Laplace's partial differential equation (Equation (1)), the Dirichlet boundary conditions imposed on the system in which the electric potential in the solution should be known and specified, while the temperature initial and Neumann boundary conditions should be known to solve the heat partial differential equation (Equation (2)). Once the temperature in a specific spatial point is obtained, Equation (4) can be solved to calculate the food new local electrical conductivity value. Accordingly, the solution of the next spatial point in the OH cell is calculated following the same previous steps.

Collinear Ohmic Heater Model
In support to the statement of Kaur N. and Singh A.K. [19] that the most commonly used OH systems in the industry are of the continuous flow type, Bozkurt H. and Icier F. [20] reported that the commercially used ohmic heating plants in pasteurization or sterilization of liquid foods (juices, milk, egg, purees, and pulps) throughout the world are of the continuous flow type.
There are two possible continuous OH configurations to be considered, the transverse configuration and the collinear configuration. The main difference between these configurations is the food product flow direction to the applied electric field where it is perpendicular in the transverse configuration and parallel in the collinear one. However, Varghese K.S. and others [2] reported that the transverse configuration has two electrical drawbacks that restricted its development, the close proximity of the live electrodes at the inlet and the outlet pipe which could cause leakage currents to earth through the processed food material, and phase-to-neutral current density in the direction of the food flow at the edges of the electrodes which could cause localized overheating and electrodes erosion due to its non-uniformity. Thus, and unlike the collinear configuration, the applications of the transverse configuration have been restricted to a specific type of liquid food.
In contrast, the collinear configuration which mainly consists of two parts, the electrode housing and the spacer tube (excellent electrical isolator). Its surface area that is in contact with the food is small such that it reduces the possibility of electrodes erosion and food overheating. Furthermore, its operational voltage can reach up to (4500 V) which is along with the collinear configuration guarantee high heating power that makes it attractive to the industrial applications It has been found that modelling of the Ohmic Heater process of the collinear configuration is of a great benefit in terms of analyzing and identifying the OH process variables and dynamic. Although the majority of the literature discusses the OH process in the batch (static) configuration, Pesso T. and Piva S. [21] provided an interesting study of a continuous Ohmic Heater of the cylindrical collinear configuration used in the sterilization of an Apricot puree. They analyzed the thermos-fluid behavior of the heater in laminar flow, calculated the resulted thermal field using an analytical solution, and studied the influence of the puree mass flow rate and the OH cell degree of insulation on the process. Based on the study by Pesso T. and Piva S. [21], a collinear Ohmic Heater mathematical model is introduced below using the numerical solution approach provided by MATLAB/Simulink.

System Geometry and Assumptions
To be modelled Ohmic Heater geometry is illustrated in Figure 1. The heater consists of a circular electrically insulating glass pipe (spacer tube) and two circular stainless steel electrodes placed at the pipe ends (inlet and outlet). The heating section (heater length) is 2.5 m, the inner diameter of the glass pipe is 50 mm, the outer diameter is 60 mm, and the applied voltage across the electrodes is 2304 V.
The working fluid is Apricot puree which has the same thermo-physical properties as water (see Table 1), its electrical conductivity (σ o ) is 0.998 S/m, and its temperature coefficient (k o ) is 0.015 1/˚C. The thermal coupling at the external surface of the pipe is considered in the modelling and the ambient temperature (T a ) is taken as 20˚C.
The working fluid is assumed to be a homogeneous one, with a 40˚C uniform inlet temperature (T in ), the current density over any cross-section of the pipe is assumed to be constant, the fluid viscous dissipation and axial heat conduction in the wall of the glass pipe are neglected, the glass pipe wall is assumed to be electrically insulated, the heating process along the heater horizontal axis is assumed to be symmetrical, and all thermos-physical properties of the fluid

System Governing Equations
Since the OH spacer tube is assumed to be electrically insulated, and the current density over any cross-section of the heater tube is constant, Equation (1) is simplified to the following form (Equation (5)) [21].
The electric potential V is calculated according to the location along the heater horizontal axis X (heater length), which is expected as the heater length is much bigger than its diameter ( L R  ), hence, the electric potential values in its Y and Z axis are constants. The system Dirichlet boundary conditions are as follows.
Equation (6a) states that the electric potential at the inlet (x = 0) of the heater is equal to the half of the applied voltage across the electrodes with a positive polarity, while Equation (6b) defines the electric potential value at the outlet (x = L) to be equal to the half of the applied voltage with a negative polarity.
Since the puree electrical conductivity is a function of its temperature and not of the heater horizontal axis X, Equation (5) can be further simplified considering the electrical conductivity as a constant to be in the following form (Equation (7)) To solve Equation (7), a MATLAB code has been developed using diff and dsolve functions.
After the electric potential (V) across the heater horizontal axis is calculated, the electric field values can be obtained, a MATLAB code has been developed to obtain the electric field values across the heater x-axis using Forward Finite Difference Method (Equation (8)) [22].
Applying the aforementioned solution, it has been found that the electric field is constant across all points of the heater x-axis with a value equal to (−921.6 V/m) (see Figure 2). The obtained solution agrees with Pesso T. and Piva S. [21] that the electrical potential in the fluid is constant with a value equal to (E = V/L).
With regard to the heater thermal field governing equation and since the working fluid is assumed to be homogenous such that it flows in the laminar pattern, the heater thermal behavior is governed by (Equation (9)).  where W is the fluid mass flow rate (kg/hour), r (m) is the heater tube radius.
With regard to the heater thermal field governing equation and since the working fluid is assumed to be homogenous such that it flows in the laminar pattern, the heater thermal behavior is governed by (Equation (9)).
where W is the fluid mass flow rate (kg/hour), r (m) is the heater tube radius.
The system Initial and Neumann boundary conditions are as described in Equations (10a), (10b) and (10c).
where H is the heat transfer coefficient (W/m 2 •℃).
Equation (9) is the same as Equation (2)  which states that the velocity of the liquid that in immediate contact with the solid (pipe wall) is the same as the velocity of the solid, which means the fluid velocity at (2r = D) is equal to zero as it can be found in the above term. Since the velocity of the fluid closer to the pipe surface is lesser than that in the center of the pipe, the fluid surface is heated more than its center which may affect the resulted heat uniformity. Equation (10b) specifies the heat flux at the heater center to be equal to zero, while Equation (10c) defines the heat conduction at the pipe surface with the ambient temperature T a . To solve Equation (9), a MATLAB function has been developed, the function uses MATLAB's Partial Differential Equation Toolbox which uses finite element analysis (method) to obtain the equation numerical solution. Since the fluid is flowing within the heater cell from the inlet to the outlet with a specific velocity determined by the fluid mass flow rate, the simulation time is equal to the fluid residence time within the heater which is calculated based on the fluid velocity (mass flow rate).
Createpde is a MATLAB function that returns the OH thermal analysis model and since the fluid temperature would increase as it flows through the heater, the analysis mode has been specified as a transient one. Generate Mesh function creates the Ohmic Heater cell mesh based on the pre-specified cell geometry, and in order to analyze the temperature along the heater diameter a max step size between any two nodes in the generated mesh has been specified to be equal to the half of the heater diameter (D/2). After the OH cell mesh is created, the fluid thermo-physical properties are specified and Equation (9) is solved according to the system initial and boundary conditions using developed MATLAB code.

MATLAB Functions and Simulink Models
Based on the aforementioned obtained and designed solutions, MATLAB functions and Simulink models have been developed as follows.

MATLAB Functions
A MATLAB function named Ohmic has been developed, the function models the continuous collinear ohmic heating process using the aforementioned solutions.
The function takes 14 process factors as inputs illustrated in Table 2.  The function defines all of its inputs as global variables to be used within other called functions. In accordance to the iterative solution in [21] the function assumes an ideal axial temperature distribution across the heater length (L) and calculates the fluid electrical conductivity according to each point temperature.
The function solves the Ohmic Heater electric field equation using MATLAB codes described earlier and solves the thermal field equation using the aforementioned solution function.
After the thermal field solution is obtained, the radial temperature distribution at the centre, halfway between the centre and the surface, and the surface are extracted using the MATLAB commands. The working fluid can be considered as a pure electrical resistance, therefore, the consumed electric power can be calculated using voltage times current. The calculation of the electrical current can be based on Ohm's law which gives the relation that current density is directly proportional to the electric field.

Simulink Model
A Collinear Ohmic Heater Simulink model has been developed using MATLAB S-function (system function) which provides a powerful mechanism for extending the capabilities of the Simulink environment, S-function has dynamically linked subroutines that the MATLAB execution engine can automatically load and execute. Basically, it provides a way to build a user-defined Simulink block that executes integrated MATLAB functions. This user-defined Simulink block is shown in Figure 3 and Fluid Velocity Sub-System model is presented as an example in Figure 4.
Ohmic Heater is the Simulink block name, setup is a function that sets the block pre-specified parameters, the block has 12 inputs (initial temperature, ambient temperature, voltage, cell length, cell diameter, fluid velocity, fluid density, thermal conductivity, specific heat, heat convection coefficient and electrical conductivity), 6 outputs (inlet temperature, current temperature, internal heat generation, electrical conductivity, electrical power and fluid position), the block inputs and outputs ports are set to be dynamic and defined to be of the Real data type, the block sampling time is specified to be continuous, and its simulation state is set to the default value.

Model Validation
The third stage of building the OH mathematical model is the validation (testing) stage. The obtained MATLAB/Simulink model is validated against the analytical study results in [21], as mentioned earlier the obtained analytical solution model outlet temperature is examined as well as the fluid mass flow rate and the cell degree of insulation influences on the process outlet temperature, axial temperature distribution, and radial temperature distribution.
In the first case study, the external convective heat transfer coefficient is estimated as 5 (W/m 2 •˚C) for natural convection and the external radiative heat Open Journal of Modelling and Simulation  transfer coefficient has been taken as 5 (W/m 2 •˚C), therefore, the heat transfer coefficient (convection and radiation) in total is 10 (W/m 2 •˚C), the case study process factors are demonstrated in Table 3 [21]. Using the above process factors values in the obtained MATLAB/Simulink model, the following outlet temperature results are obtained in Table 4.
The results of developed MATLAB/Simulink model are not identical to the analytical solution provided by Pesso T. and Piva S. [21] which is expected as the numerical approach provides an approximate solution to the problem, but it showed good outlet temperature estimations with an acceptable error range. Note that the obtained results can be further enhanced by decreasing the step size of the cell generated mesh, for instance, and for the same case study with a (D/8) mesh step size, the MATLAB/Simulink model results were improved (99.97˚C) however, the computation time increased significantly.
In the second case study, Pesso T. and Piva S. [21] examined the fluid mass flow rate influence on the process outlet temperature, using the same process factors (Table 3) but the convection heat transfer coefficient. Its value is calculated according to Equation (11).
where K is the thermal conductivity (W⁄m•˚C), D i is the heater diameter (m), and H i is the convection heat transfer coefficient (W/m 2 •˚C).
The study uses R w (thermal resistance) with a value of 1.1 which is the equivalent to 11.233 W/m 2 •˚C, and since the external radiative heat transfer is 5 W/m 2 •˚C the total heat transfer coefficient value is 16.233 W/m 2 •˚C. Table 5 illustrates the obtained model's results in comparison with the study results in [21].   The model yielded good results in comparison with analytical solutions in [21] with an acceptable error range which could be decreased with a smaller mesh step size, however, and for this study purposes, the generated solutions are acceptable. In the third case study, Pesso T. and Piva S. [21] examined the thermal insulation influence on the process outlet temperature, using the same aforementioned process factors in Table 3, the equivalent convection heat transfer coefficient values to the thermal resistances provided in the study are calculated based on Equation (11) plus the value of the external radiative heat transfer. Table 6 demonstrates the obtained model's results in comparison to the study outcomes.
The model showed the same variable influence as the analytical solution, yet its outlet temperature errors in comparison to results reported in [21] are slightly larger than in the previously discussed case studies. Figure 5 illustrates MATLAB/Simulink model axial temperature distribution, which shows its non-linearity as a function of the heater length.
The model showed that the temperature is non-linearly distributed along the heater horizontal axis. This result agrees with Pesso T. and Piva S. study in [21] that the bulk temperature distribution is non-linear. Figure 6 shows the analytical solution temperature distribution in [21].
In comparison with the obtained model's results, the axial temperature   distribution is non-linear, yet the model's outputs slopes are slightly different than the analytical model slope as seen in the above figures. The radial temperature profile produced by the obtained model has the same pattern as the analytical solution in [21] (Figure 7 & Figure 8).
The developed model has the same radial temperature pattern as the analytical model with some temperature differences, besides that the model's estimation is approximate and not exact, the radial temperature profile has been calculated for three points only (r= 0, r= 0.0125, and r = 0.025) which yielded some estimations differences.
According to the aforementioned case studies the MATLAB/Simulink model temperature estimation accuracy range is 97.8-99.6% (Tables 4-6) therefore, and in accordance with the study purposes in terms of process analysis and control systems design, the models are found to be valid and acceptable for further use.

Conclusions
Mathematical models are invaluable tools that describe systems' behaviour in the real world and through which systems are investigated and analysed. The Ohmic Heating process consists of two sub-processes, the electric process, and the thermal process, both processes governing equations are of the partial differential types which require boundary and initial conditions to be solved.
The analytical and the numerical solutions are methods used to obtain OH governing equations solutions, the analytical solution is highly accurate, while the numerical solution provides approximate results, however, the numerical solution can be achieved using computer-aided software. The Collinear Ohmic Heater configuration is the most commonly used configuration in the industry.
The mathematical model of its heating process has been developed using MATLAB/Simulink. The obtained model showed good output results in terms of outlet temperature, axial temperature distribution and radial temperature distribution with an accuracy up to 99.6% in comparison to the analytical solution.