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  .
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 . 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  .
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  . Goullieux and Pain  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.  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 19th 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 . 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 . Sarang and others  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  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 , as the Ohmic Heating process has close to 100% energy transfer efficiency, while Microwave heating has 65% efficiency at its best .
Kamonpatana and others  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  .
Ohmic Heating is one of the emerging thermal technologies that could be better developed, understood and validated through its mathematical model . 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 .
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  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  modelled a 3-dimensional OH for high frequencies based on Maxwell’s equations of the electric field, while De Alwis and Fryer  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.
2. Mathematical Modelling of Ohmic Heating Processes
2.1. 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 .
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 .
where ρ is the food density (kg/m3), Cp 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/m3) 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) To, ko 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, the system boundary and initial conditions should be specified. Mainly, there are three types of boundary conditions; Dirichlet boundary condition, Neumann boundary condition, and Mixed boundary condition. Dirichlet boundary condition specifies the value of the function at the boundary of the system geometry, Neumann boundary condition defines the value of the function derivative at the boundary of the system geometry, while the Mixed boundary condition forms a combination of both .
After the system boundary and initial conditions are specified, the system PDEs could be solved. Generally, there are two main solutions through which partial differential equations are solved, the analytical solution and the numerical solution. In terms of the solution accuracy and the understanding of the presented problem the analytical solution is usually preferred, however, when models contain non-linear components analytical results are difficult to obtain, thus, a numerical solution is often used.
The numerical solution mimics the model process and provides an approximate solution to the presented problem, there are multiple methods through 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.
2.2. Collinear Ohmic Heater Model
In support to the statement of Kaur N. and Singh A.K.  that the most commonly used OH systems in the industry are of the continuous flow type, Bozkurt H. and Icier F.  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  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.  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. , a collinear Ohmic Heater mathematical model is introduced below using the numerical solution approach provided by MATLAB/Simulink.
2.2.1. 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 (ko) is 0.015 1/˚C. The thermal coupling at the external surface of the pipe is considered in the modelling and the ambient temperature (Ta) is taken as 20˚C.
The working fluid is assumed to be a homogeneous one, with a 40˚C uniform inlet temperature (Tin), 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
Figure 1. Collinear ohmic heater.
Table 1. Apricot puree thermo-physical properties.
except for the electrical conductivity are assumed to be constant (not changing with the temperature).
2.2.2. 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)) .
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 ( ), 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)) .
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.  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)).
Figure 2. Electric field and voltage potentials.
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/m2∙℃).
Equation (9) is the same as Equation (2) but in the cylindrical coordinates, the term expresses the flowing fluid no-slip boundary condition
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 Ta.
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.
3. MATLAB Functions and Simulink Models
Based on the aforementioned obtained and designed solutions, MATLAB functions and Simulink models have been developed as follows.
3.1. 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.
Table 2. Ohmic function inputs.
The function defines all of its inputs as global variables to be used within other called functions. In accordance to the iterative solution in  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.
3.2. 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.
OhmicHeater 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.
4. 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 , 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/m2∙˚C) for natural convection and the external radiative heat
Figure 3. Ohmic heater simulink block.
Figure 4. Fluid velocity sub-system.
transfer coefficient has been taken as 5 (W/m2∙˚C), therefore, the heat transfer coefficient (convection and radiation) in total is 10 (W/m2∙˚C), the case study process factors are demonstrated in Table 3 .
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.  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.  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), Di is the heater diameter (m), and Hi is the convection heat transfer coefficient (W/m2∙˚C).
The study uses Rw (thermal resistance) with a value of 1.1 which is the equivalent to 11.233 W/m2∙˚C, and since the external radiative heat transfer is 5 W/m2∙˚C the total heat transfer coefficient value is 16.233 W/m2∙˚C. Table 5 illustrates the obtained model's results in comparison with the study results in .
Table 3. First experiment factors values.
Table 4. OH model results.
Table 5. OH model results.
Table 6. OH model results.
The model yielded good results in comparison with analytical solutions in  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.  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  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  that the bulk temperature distribution is non-linear. Figure 6 shows the analytical solution temperature distribution in .
In comparison with the obtained model’s results, the axial temperature
Figure 5. MATLAB/Simulink model axial temperature distribution.
Figure 6. Analytical model axial temperature distribution .
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  (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.
Figure 7. Matlab/simulink model radial temperature distribution.
Figure 8. Analytical model radial temperature distribution .
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. MATLAB/Simulink model is found to be valid for further studies and useful in control applications where real-time simulation of the process is required.
The authors would like to thank H2020 ERA-NET and DEFRA for providing funding for this research. The authors would also like to thank Sheffield Hallam University for providing the APC.
 Goullieux, A. and Pain, J. (2014) Ohmic Heating. In: Sun, D.W., Ed., Emerging Technologies for Food Processing, Academic Press, Amsterdam, 399-420.
 Varghese, K.S., Pandey, M.C., Radhakrishna, K. and Bawa, A.S. (2014) Technology, Applications and Modelling of Ohmic Heating: A Review. Journal of Food Science and Technology, 51, 2304-2317.
 Kamonpatana, P., Mohamed, H., Shynkaryk, M., Heskitt, B., Yousef, A. and Sastry, S. (2013) Mathematical Modeling and Microbiological Verification of Ohmic Heating of a Solid-Liquid Mixture in a Continuous Flow Ohmic Heater System with Electric Field Perpendicular to Flow. Journal of Food Engineering, 118, 312-325.
 Simpson, R., Nuñez, H., Jaques, A., Ramírez, C., Quiroz, N., Moreno, J. and Sastry, S. (2018) Application of a Moderate Electric Field for the Potential Acceleration of the Salting Process of Atlantic Salmon (Salmo salar). Journal of Food Process Engineering, 41, e12846.
 Gotz, A., Wani, A.A., Langowski, H.C. and Wunderlich, J. (2014) Food Technologies: Aseptic Packaging. In: Yasmine, M., Gerald, M. and Ewen, T., Eds., Encyclopedia of Food Safety, Academic Press, Cambridge, 124-134.
 Palaniappan, S. and Sastry, S.K. (1991) Electrical Conductivities of Selected Solid Foods during Ohmic Heating. Journal of Food Process Engineering, 14, 221-236.
 Jaeger, H., Roth, A., Toepfl, S., Holzhauser, T., Engel, K.H., Knorr, D., Vogel, R.F., Bandick, N., Kulling, S., Heinz, V. and Steinberg, P. (2016) Opinion on the Use of Ohmic Heating for the Treatment of Foods. Trends in Food Science & Technology, 55, 84-97.
 Sarang, S., Sastry, S.K. and Knipe, L. (2008) Electrical Conductivity of Fruits and Meats during Ohmic Heating. Journal of Food Engineering, 87, 351-356.
 Cho, W.I., Kim, E.-J., Hwang, H.J., Cha, Y.H., Cheon, H.S., Cho, J.B. and Chung, M.S. (2017) Continuous Ohmic Heating System for the Pasteurization of Fermented Red Pepper Paste. Innovative Food Science & Emerging Technologies, 42, 190-196.
 Marra, F., Zell, M., Lyng, J.G., Morgan, D.J. and Cronin, D.A. (2009) Analysis of Heat Transfer during Ohmic Processing of a Solid Food. Journal of Food Engineering, 1, 56-63.
 Shim, J.Y., Lee, S.H. and Jun, S.J. (2010) Modeling of Ohmic Heating Patterns of Multiphase Food Products Using Computational Fluid Dynamics Codes. Journal of Food Engineering, 99, 136-141.
 Guo, W., Llavea, Y., Jin, Y., Fukuoka, M. and Sakai, N. (2017) Mathematical Model of Ohmic Heating of Two-Component Foods with Non-Uniform Electric Properties at High Frequencies. Innovative Food Science and Emerging Technologies, 39, 63-78.
 De Alwis, A.A.P. and Fryer, P.J. (1990) A Finite Element Analysis of Heat Generation and Transfer during Ohmic Heating of Foods. Chemical Engineering Science, 45, 1547-1559.
 Jun, S. and Sastry, S. (2005) Modeling and Optimization of Ohmic Heating of Foods Inside a Flexible Package. Journal of Food Process Engineering, 28, 417-436.
 Kumar, M. and Mishra, G. (2011) An Introduction to Numerical Methods for the Solutions of Partial Differential Equations. Applied Mathematics, 2, 1327-1338.