- Email: [email protected]
Freezing times using time derivative of temperature on surface of foods
Accepted Manuscript
Freezing times using time derivative of temperature on surface of foods S.R. Ferreira , L.O.A. Rojas PII: DOI: Reference:
S0140-7007(18)30442-0 https://doi.org/10.1016/j.ijrefrig.2018.11.009 JIJR 4168
To appear in:
International Journal of Refrigeration
Received date: Revised date: Accepted date:
2 June 2018 26 October 2018 7 November 2018
Please cite this article as: S.R. Ferreira , L.O.A. Rojas , Freezing times using time derivative of temperature on surface of foods, International Journal of Refrigeration (2018), doi: https://doi.org/10.1016/j.ijrefrig.2018.11.009
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Highlights A numerical model was developed to predict the food freezing process. (FDM) was employed to calculate the temperature distribution in a food T(x,t). This paper tests the use of dT(i=n)/dt in x=L, and dT(i=1)/dt in x=0. Simulations obtaining (tcalc) were compared with 227 experimental data (texper).
AC
CE
PT
ED
M
AN US
CR IP T
Good agreement was obtained between predicted and experimental values.
ACCEPTED MANUSCRIPT Freezing times using time derivative of temperature on surface of foods S.R. Ferreira a*, L.O.A. Rojas b a
Department of Chemical Engineering, UFRN - University Federal of the Rio Grande
do Norte, Campus Universitario - 59078-970 - Natal - RN - Brazil. Tel. 55-8432153753 b
Department of Chemical Engineering, UFPb – Federal University of the Paraiba -
*
[email protected]
Abstract
CR IP T
Joao Pessoa - Pb - Brazil.
A numerical model is developed to predict the freezing process of foods or analogues.
AN US
The finite difference method (FDM) is employed to calculate the temperature
distribution in a food T(x,t). The aim of this research study is to test the use of the time derivative of temperature on the surface, dT(i=n)/dt in x=L, and at the center of food, dT(i=1)/dt in x = 0. This methodology is employed to predict the freezing time (tcalc),
M
using the properties of tylose gel (k, ρ, and Cp), which is a food analog. Simulations were performed using the model, obtaining (tcalc), and compared with 227 experimental data points (texper). In summary, the main conclusions of this study are as follows. It is
ED
important to use the derivative dT(n)/dt at x = L in model (2) in relation to the use of dT(1)/dt at x = 0 in the same model. The evaluation of thermal conductivity using
PT
temperature T = Tn-2, k(Tn-2) at x = L with the proposed model and properties of tylose gel to calculate (tcalc) is adequate. Simulations using k(Tn-2) at x = L with the proposed
CE
model were performed, obtaining (tcalc), and compared with 227 experimental data points (texper), exhibiting good agreement. The following parameters were obtained for
AC
227 experimental data points: minimum error Emin = -6.97 %, overall mean error Emean = 0.06 %, maximum error Emax = 8.08 %, and standard deviation n-1 = 2.19 %. Keywords: time derivative of temperature; convective boundary condition; freezing time; finite difference method.
Nomenclature
Bi
Biot number [--]
ACCEPTED MANUSCRIPT Bif
Biot number in the fully frozen state [--]
Bim
mean value of Biot numbers of the sample in the fully frozen state [--]
Bimax maximum value of Biot number in the fully frozen state [--] minimum value of Biot number in the fully frozen state [--]
Cp
effective specific heat [J kg-1 °C-1]
Cp0
unfrozen effective specific heat [J kg-1 °C-1]
Cp1
freezing zone effective specific heat [J kg-1 °C-1]
dx
step size for spatial variable [m]
E
error in calculation or experimental percentage error [%]
h
convective heat transfer coefficient [W m-2 °C-1]
H
specific enthalpy [J kg-1]
i
ith node [--]
k
thermal conductivity [W m-1 °C-1]
kf
thermal conductivity of food in fully frozen state [W m-1 °C-1]
k(i)
thermal conductivity calculated at temperature of ith node T(i) [W m-1 °C-1]
k0
unfrozen thermal conductivity [W m-1 °C-1]
k1
freezing zone thermal conductivity [W m-1 °C-1]
L
slab half-thickness [m]
n
number of nodal points or on surface [--]
np
number of divisions in time [--]
R
radius of infinite cylinder or sphere, or object half-thickness (L) [m]
t
time [s]
tcalc
calculated freezing time [s]
texper
experimental freezing time [s]
AN US
M
ED
PT
integration time [s] temperature [°C]
AC
T
CE
tend
CR IP T
Bimin
Ta
ambient or cooling medium temperature [°C]
Tc
center temperature or desired final (T) at center of food [°C]
Tf
initial freezing temperature [°C]
Tm
average temperature at surface in i = n for k(Tm) calculation [°C]
Tn
temperature at surface with convection [°C]
T0
initial temperature [°C]
x
distance measured along the x-axis [m]
Xmax
maximum deviation [%]
ACCEPTED MANUSCRIPT Xmean mean deviation [%] Xmin
minimum deviation [%]
Y
temperature in computational program [°C]
Ya
cooling medium temperature in computational program [°C]
Ym
average temperature at surface in i = n for k(Tm) calculation [°C]
Greek symbols
distance increment in (x) dimension [m]
ρ
density [kg m-3]
ρ0
unfrozen density [kg m-3]
ρ1
freezing zone density [kg m-3]
σn-1
standard deviation [%] or [--]
M
Subscripts
AN US
∆x
ambient
c
center
calc
calculated value
end
end time of calculation in computational program
ED
a
PT
exper experimental value freezing
i
ith
m
mean or average maximum
AC
max
CE
f
min
minimum
n
number of nodes in (x) dimension
0
initial or unfrozen
1
in freezing zone
1. Introduction
CR IP T
Yprime time derivative T/t in computational program [°C s-1]
ACCEPTED MANUSCRIPT Freezing times (tcalc) are estimated to predict the required residence times for foods in freezers. As foods freeze, a frozen layer forms next to the outer surface and moves progressively inward. Ice grows, and heat is removed, mainly at, or near the inwardmoving freezing interface (Schwartzberg et al., 2007, p. 70). Commercial finite element method (FEM) software often provide graphical user interfaces, are user-friendly, can be used without understanding FEM principles, and have been used to treat complex and noncomplex freezing and thawing problems
CR IP T
(Schwartzberg et al., 2007, p. 89).
Several numerical methods have been developed to solve partial differential
equations (Northrop et al., 2013). For parabolic equations such as the heat equation,
there exist several numerical methods to find a solution (Dehghan, 2006; Northrop et
al., 2013). The method of lines (MOL) is one such efficient routine in which the space is
AN US
discretized using any of a number of techniques such as the finite difference method (FDM), finite element method (FEM), finite volume method (FVM), or collocation method (CM) (Cutlip and Shacham, 1999; Dehghan, 2006; Northrop et al., 2013;
Sadiku and Obiozor, 2000; Schiesser, 1991; Schiesser and Silebi, 2009). Ferreira et al.
M
(2016) and Ferreira (2017) used the MOL to predict (tcalc), based on the effective specific heat method.
Generally, commercial software based on the FDM are unable to run and/or converge
ED
to accurate temperature predictions when attempting to simulate a phase change process, owing to abrupt changes in the thermophysical properties, specifically when
PT
using the apparent specific heat method (Dima et al., 2014). The surface of an object is the most difficult to model in the freezing process when
CE
there is a convective boundary condition (Albasiny, 1956; Bonacina and Comini, 1973, p. 586-587; Cleland, 1977, p. 18-19; Cleland and Earle, 1977b, p. 1030-1032). For this
AC
reason, the focus of this approach is to determine the influence of the time derivatives of temperature, dT(i=n)/dt and dT(i=0)/dt. This approach allows the calculation of the heat transfer convective coefficient (h) and an accurate prediction of the (tcalc). Based on what has been cited, the objectives of this research are as follows: To develop a mathematical model, called model (2), for simulating the freezing process in a slab, and calculate the (tcalc). To evaluate the influence of the use of dT(i=n)/dt and dT(i=0)/dt. To compare the predictions of the proposed model (tcalc) to published experimental data (texper).
ACCEPTED MANUSCRIPT Experimental freezing times were selected from published data (Bazan and Mascheroni, 1984; Cleland and Earle, 1977a; Creed and James, 1983, 1985; Hense and Kieckbusch, 1991, 1985; Hung and Thompson, 1983; Lescano, 1973; Mascheroni and Calvelo, 1982; Michelis and Calvelo, 1983; Pham and Willix, 1990).
2. Literature review
CR IP T
In this section, some published numerical methods that are used to predict (tcalc), and several challenges of the phase change process modeling are briefly reviewed.
2.1. Some numerical methods for freezing-time prediction
AN US
A survey of the literature since 1984 shows that this systematic approach for assessing prediction methods has not been widely used (Cleland et al., 1994). This survey also shows that no method (including numerical methods) predicts all experimental data well, indicating that some data that are widely used for testing are of lower quality than
M
that of others. Furthermore, it may be possible to ensure that the data collected in the future is more accurate (Cleland et al., 1994). The choice of technique depends on the objectives of the modelers and the technical
ED
means at their disposal (Pham, 2008, 2014). Numerical methods can, in principle, provide exact or near-exact predictions for almost any situation, although in practice,
PT
their accuracy is limited by an inadequate knowledge of the problem parameters, for example, the product properties, geometry, and flow characteristics (Pham, 2008, 2014).
CE
Schwartzberg et al. (2007) discussed the apparent specific heat formulation and the enthalpy step method for freezing and thawing calculations, including the numerical
AC
solution of the partial differential equations (PDEs) describing the heat transfer during these processes. One such method, based on the use of specific heat (Cp), calculates the T(x,t) using
explicit methods, such as Euler‟s forward difference method. The following information provided by Schwartzberg et al. (2007) is relevant for this method of calculating the (tcalc): The (tcalc) and T(x,t) for food freezing, obtained using these methods, are presented.
ACCEPTED MANUSCRIPT The computation becomes unstable if the time step (t) is too large. Stability is provided by using the smallest (t) obtained from the equations presented in the literature (Mannapperuma and Singh, 1988, 1989; Schwartzberg et al., 2007). When the equations for obtaining T(x,t) were solved for typical freezing conditions by machine computation, 700–2,000 time steps were used for slabs, 370–1,000 steps for infinite cylinders, and 130–1,500 steps for spheres. The number of steps needed is inversely proportional to (L2) and proportional to
CR IP T
n2(n+1), but also depends on Bif = 2hL/kf, (T0), (Ta), and (Tc); (n) is the number of
nodal points for FDM calculations, (Bif) is Biot number in the fully frozen state, (h) the convective heat transfer coefficient, (L) the object‟s half thickness or radius (R), (kf) the thermal conductivity of the food in the fully frozen state, (T0) the initial
AN US
temperature, (Ta) the temperature of the cooling, and (Tc) the desired final (T) at the center of the food. The typical Biot number is defined as Bi = hL/k (Luikov, 1968; Singh and Heldman, 2009); which is different from that employed by Schwartzberg et al. (2007) and Pham and Willix (1990).
One can use a much larger (t) with both the Lee and Crank–Nicolson methods than
M
that for the explicit method, but if (t) is too large, „jumping‟ and (T) oscillations
ED
occur as (T) passes through (Tf), where (Tf) is the initial freezing point.
2.2. Challenges of the phase change process modeling
PT
The problem of predicting the time for a material to freeze under a given set of conditions is classed physically as one of heat conduction with change of phase
CE
(Cleland, p. 18, 1977).
The change of phase can be taken into account in two different ways. It can be
AC
assumed that all the latent heat is released at a freezing temperature at which a step change in thermophysical properties occurs (Carslaw and Jaeger, p. 282, 1959; Cleland, p. 18, 1977). Alternatively, the release of latent heat over a range of temperatures can be taken into account by changing specific heat Cp(T) and thermal conductivity k(T). The first has been used extensively, even though variations in thermal properties mean that a freezing foodstuff will fit the solution only approximately. The second method introduces a problem in defining the completion of freezing, which is deemed complete when the phase change front reaches the thermodynamic centre of the material in the
ACCEPTED MANUSCRIPT first method, and the presence of a mushy zone precludes the use of such a definition (Cleland, p. 19, 1977). There are several challenges for modelers in the food freezing process, such as trying to understand the phase change, mathematical modeling of the process and trying to solve the numerical problem to calculate the freezing time. The freezing of foodstuffs frequently involves boundary conditions of the third kind, which have proved to be difficult to handle and subject to error, mainly because of the
CR IP T
requirement to calculate surface temperatures implicitly (Cleland and Earle, p. 1029, 1977b).
The focus of this section is the boundary condition on the surface of a slab, because it is the most difficult point to be modeled, that is, to analyze the use of derivative dT(i=n)/dt. The use of derivative dT(i=1)/dt in x = 0 is also tested.
AN US
In design situations the surface temperature/time profile is rarely known. However
the bulk average ambient temperature (Ta) is often known, for example, air temperature in a blast freezer. There are also well known means of arriving at a surface heat-transfer coefficient (h), which considers all resistance for heat leaving the surface of the freezing
M
material and entering the cooling medium (Cleland and Earle, p. 1030, 1977b). To facilitate analysis, a slab is considered. The half-thickness of the slab (L) in the range of 0 ≤ x ≤ +L is divided into (n-1) sections with (n) node points, numbered from
ED
(1) to (n); that is, setting i = 1 in the center of the slab when x = 0, and i = n on the slab surface when x = L. The space grid is shown in Fig. 1.
PT
In these cases, the surface requires special consideration. To solve some of the difficulties and improve freezing time accuracy, a number of researchers proposed
CE
changes in the mathematical modeling on the surface of a slab. Some of these are described below:
AC
Bonacina and Comini (1973, p. 586) suggest that when a fictitious point i = (n+1) is used, then T(n+1) = Ta, and pseudo-thermal conductivity is defined, k(n+1/2) = hx, so that a three level scheme can be used. The authors found the “Lees” scheme to be superior to others, due to smaller truncation errors, and better representation of thermal properties (Cleland and Earle, 1977b, p. 1030). It is important to underscore, however, that the three-time level method remains unconditionally stable, even with boundary conditions of the second and third kind (Bonacina and Comini, 1973, p. 587).
ACCEPTED MANUSCRIPT In another approach, the finite difference boundary condition was derived by Cleland and Earle (1977b, p. 1032) from a heat balance over the surface space increment which extends (0.5x) from the surface. While the internal nodes have the specific heat (Cp) of a full (x) associated with them, the surface node has the (Cp) of only a half space increment. The freezing times calculated from this approach were compared to their experimental data (Cleland, 1977, p. 103).
CR IP T
Albasiny (1956) used the enthalpy method and a fictitious point i = (n+1) with T(n+1) = Ta, obtaining a time derivative of enthalpy on the surface of the slab,
dH(i=n)/dt in x=L. The set of resulting equations was solved using the Runge-Kutta method, obtaining enthalpy (H) and the corresponding temperature (T) versus time (t). According to Cleland (1977, p.97) the solution obtained by Albasiny (1956) is
due to the large truncation error.
3. Materials and methods
AN US
not as accurate as the scheme of Lees (1966) or the Crank-Nicolson (1947) method
M
The method is described in five sections. The freezing process model is described in Section 3.1, based on the heat conduction equation for a slab with the denominated
ED
model (2). This model is discretized using an FDM, as shown in Section 3.2, and solved using a computational program by the MOL, as shown in Section 3.3. The selected experimental (tcalc) dataset from the literature is presented in Section
PT
3.4. Finally, the selection of the physical properties that are employed as parameters in the numerical program when implementing model (2) is discussed in Section 3.5. The
CE
results obtained are discussed in the Results and discussion section.
AC
3.1. Unsteady-state freezing for symmetric heat conduction in slabs
For illustrative purposes, the freezing process in slab-shaped foods is analyzed. Analogous equations can be written for spheres and infinite cylinders. A slab located within the range of -L ≤ x ≤ +L is considered. When the heat transfer is symmetric, it is convenient to use the slab center as the origin for the coordinates. Because of symmetry, only the region between x = 0 and x = L is treated. The half-
ACCEPTED MANUSCRIPT thickness of the slab (L) considered is divided into (n-1) sections with (n) node points, indexed from (1) to (n) from the center to the slab border. A heat balance over a small section of a slab of thermal conductivity k(T), assuming that the material within the ith element has a specific heat (Cp) and density (ρ), can be represented by (Cleland, 1977; Cleland and Earle, 1984)
()
()
( ( )
⁄
)
⁄
At the limit when ∆x 0, Eq. (1) gives
()
()
( ( )
)
⁄
)
⁄
(1)
(2)
AN US
()
( ( )
CR IP T
()
The first term of Eq. (1) or (2), that is, the accumulation of heat in the slab, can be calculated at a temperature (Ti) with Cp(Ti) and (Ti). The last term, representing the heat conduction in the slab, can be evaluated at positions i+1/2 and i−1/2. There are
M
various options for choosing the (T), for which the properties (k), (), and (Cp) are calculated. For example, the second term of Eq. (2) can be calculated at the “centered”
ED
position (i), k(i) = k(Ti).
For a slab with thickness (2L) that is cooled from both sides at a cooling temperature
PT
(Ta), the initial and boundary conditions are
CE
(3)
AC
( )
( )
(
(4)
)
(5)
Model (2) is discussed as follows. The importance of dT(i=1)/dt and dT(i=n)/dt for the freezing-process modeling is tested. The procedure used is as follows: The slab in the range 0
x
L of Eq. (2) is discretized; this is valid from i = 1 to i =
n. However, it is usually discretized only from i = 2 to i = n-1, using the boundary
ACCEPTED MANUSCRIPT conditions at the extremes, i = 1 and i = n, to complete the energy balance for the slab at 0
x
L.
In this paper, Eq. (2) is also discretized for i = 1 and x = 0. The boundary condition of Eq. (4) is discretized and substituted into Eq. (2). It results in a discretized form for i = 1, with the first member equal to dT(i=1)/dt. Additionally, Eq. (2) is discretized for i = n and x = L. The boundary condition of Eq.
with the first member equal to dT(i=n)/dt.
CR IP T
(5) is discretized and substituted into Eq. (2). It results in a discretized form for i = n, With this procedure, all the interior and extreme points of the discretized equations will have a "local equation" as a function of dT(1)/dt, dT(2)/dt, …, dT(n-1)/dt, and dT(n)/dt.
AN US
The resulting system of equations is solved using the initial and boundary conditions. Solving the system yields the temperatures T(i) versus time (t), that is, T(1,t), T(2,t), …, T(n-1,t), and T(n,t).
M
3.2. Discretization of model (2) for slabs
After realizing the discretizations for the space, the resulting system of ordinary
ED
differential equations (ODE) is solved by the MOL. The half-thickness of the slab (L) in the range of 0 ≤ x ≤ +L is divided into (n-1)
CE
given by
PT
sections with (n) node points, numbered from (1) to (n); therefore, the space grid is
(6)
)
AC
(
Setting i = 1 in the center of the slab when x = 0, and i = n on the slab surface when x
= L, Eq. (3) can be written as ()
( )
(7)
Eq. (2) can be discretized using a central second-order approximation around point (i). The properties ρ(i), Cp(i), and k = k(i) are evaluated at T(i) for the interior points 2 ≤ i ≤ n-1, resulting in
ACCEPTED MANUSCRIPT
( )(
()
) ()
( )(
(8)
)
Salvadori (1994) presents a similar space discretization for a slab, although they performed a dual time and space discretization. However, the temperature (T) at which (k) is calculated varies from author to author; for example, Cleland (1977, p. 108) and
CR IP T
Cleland and Earle (1979, p. 962) used k(i) = (1/2)(k(i-1/2)+k(i+1/2)), while Wilson and Singh (1987, p. 151) used k(i-1/2) = k((1/2)(T(i)+T(i-1))) and k(i+1/2) = k((1/2)(T(i)+T(i+1))) etc.
(
)
(
)( ( ) ( )(
) )
AN US
Setting i = 1 in the center of the slab when x = 0, Eq. (2) can be discretized as
(9)
The properties ρ(1), Cp(1), and k = k(1) were evaluated at temperature T(1) using Eq. (9). The subsequent evaluation of the properties ρ(1), Cp(1), and k = k(2) could also be
M
tested, considering that (k) is evaluated at a mean temperature T(2) etc. Various discretization schemes can be tested for the heat flow qx = -kdT/dx for x = 0
)
(10)
PT
( )( ( )
ED
using Eq. (4) (Cutlip and Shacham, 1999, p. 396). Discretizing Eq. (4) at i = 1,
CE
Rearranging Eq. (10),
AC
(11)
Substituting Eq. (11) in Eq. (9),
(
)
( )( ( )
) ( )(
)
( )( ( )
)⁄( ) ( )( )
(12)
ACCEPTED MANUSCRIPT The energy balance is complete at x = 0 when substituting Eq. (11) in Eq. (9), which results in Eq. (12), that is, T(i=1)/t. This is due to the symmetry of the temperature profile T(x,t) at i = 1, resulting in qx = 0 at i = 1, that is, T(i=1)/t at x = 0 has less influence on the calculations of T(x,t) when compared with that of T(i=n)/t on the surface of the slab at x = L. Setting i = n on the slab surface when x = L, Eq. (2) can be backward discretized as (
)
) (
(
) (
)
( )
)
CR IP T
(
(13)
The properties (, Cp, and k) can be evaluated at a selected (T) using Eq. (13). For example, (n), Cp(n), and k(n); (n), Cp(n), and k(n-1); (n-1), Cp(n-1), and k(n-2);
AN US
(n-2), Cp(n-2), and k(n-2) etc.
Several discretizations can be tested for the heat flow dT/dx at x = L using Eq. (5) (Cutlip and Shacham, 1999, p. 396). Discretizing the boundary condition of Eq. (5) at i
(
)
( )(
)
)
(14)
ED
(
M
= n at x = L gives
)
(
)( )(
)
(15)
CE
(
PT
Rearranging Eq. (14),
AC
Substituting Eq. (15) in Eq. (13), (
)
(
)
)( )(
( (
)
)
(16)
The energy balance is complete at x = L when substituting Eq. (15) in Eq. (13), which results in Eq. (16). Generally, T(i=n)/t at x = L is more important in calculations with a convective boundary condition in comparison with T(i=1)/t at x = 0, as already mentioned. To evaluate (k) and calculate () and (Cp) to use in Eq. (16), Tm = Tn-1, Tm = Tn-2, Tm
ACCEPTED MANUSCRIPT = (1/2)(Tn-2 +Tn-1) etc. are tested, where (Tm) is the average temperature in the evaluation of k(Tm) at i = n. Eqs. (6), (7), (8), (12), and (16) can be numerically integrated to any time (t) with one adequate solver. Solving the system yields temperature T(i) versus time (t), with i=1 to i=n. The number of discretization points (n) for each experiment is found by trial and error, until the lowest percentage error E(%) is obtained between the experimental
CR IP T
(texper) and calculated values (tcalc).
3.3. Development of a computational program for slabs based on model (2)
A computational program and auxiliary subroutine were implemented in Fortran 90 to integrate the resulting system of equations, as shown in Table 1, for model (2).
AN US
The program was implemented in Intel Visual Fortran, by employing the IMSL -
Fortran Numerical Library (IMSL, 2006), using routine IVMRK. All the numerical runs were tested on a Compaq Presario CQ40 Notebook PC equipped with an Intel Core 2 Duo T6500 processor with a speed of 2.1 GHz, 4 GB of RAM, and a 64-bit operating
M
system.
ED
3.4. Selected experimental freezing time (texper) dataset from the literature Certain freezing time experimental data were selected for simulation from published
PT
data (Bazan and Mascheroni, 1984; Cleland and Earle, 1977a; Creed and James, 1983, 1985; Hense and Kieckbusch, 1991; Hung and Thompson, 1983; Lescano, 1973;
CE
Mascheroni and Calvelo, 1982; Michelis and Calvelo, 1983; Pham and Willix, 1990). Several of these experimental data are cited and analyzed in the work of Salvadori
AC
(1994).
Each freezing time value that was predicted with the proposed model (tcalc) was
compared with published experimental data (texper). The results of the calculations using model (2) are discussed in the Results and discussion section.
3.5. Process of selection of physical properties used in calculations
ACCEPTED MANUSCRIPT The set of physical properties of tylose gel (k, ρ, and Cp) employed in the freezing simulation using model (2) is called case (B) (Succar and Hayakawa, 1983; Ferreira, 2017).
4. Results and discussion
CR IP T
4.1. Simulations using experimental data from Creed and James
Table 1 illustrates the computational program and subroutine for the calculation of the slab temperature T(x,t) and freezing time (tcalc) by employing model (2).
In the 44 freezing experiments on pig liver slabs conducted by Creed and James
temperature of Tc = -7 °C is reached.
AN US
(1983; 1985), the freezing time (tcalc) is the time from the onset of cooling, until a center Table 1 presents the input parameters that are necessary for calculation with the first experiment from Creed and James (1983; 1985). The parameters for the pig liver slab are n = 7, np = 5000, tend = 10000 s, L = 0.0380 m, Y(1:n) = 10.5 °C, Ya = -21 °C, and h
M
= 361 W m-2 °C-1. The time needed to reach a value of Tc = -7 °C at the center point of the pig liver slab was calculated, resulting in tcalc = 9454.0 s = 2.63 h < tend = 10000 s. The experimental time texper = 2.57 h, resulting in a percent relative error E = 2.18 %.
ED
The run time for the program is 2 s.
The criterion employed by Pham and Willix (1990, p. 1429) was used, for which the
PT
Biot number in the frozen state Bif = h(2L)/kf is considered large if Bif ≥ 10; kf is the frozen thermal conductivity. The definition Bif = h(2L)/kf is used by numerous
CE
researchers in their studies on food freezing processes (Cleland, 1977; Pham and Willix, 1990; Schwartzberg et al., 2007). Pham and Willix (1990) used a kf = 1.54 W m-1 °C-1
AC
for tylose gel. Using this value, a Bif = 17.82 ≥ 10 is obtained for the original Experiment 1 of Creed and James (1983; 1985). The values for all the other experiments are Bif ≥ 10 (Creed and James, 1983; 1985). Table 2 presents a summary of the comparisons of the experimental and predicted
freezing times for the selected case (B), using model (2). These are some of the tests performed using model (2): Analysis of the importance of dT(1)/dt at x = 0 and dT(n)/dt at x = L in the computer program of model (2).
ACCEPTED MANUSCRIPT Comparison of the use of k(Tn-2) for the evaluation of k at x = L in relation to the use of k((1/2)(Tn-2+Tn-1)) in model (2). The relevant aspects of the prediction of (tcalc) using the data from Creed and James (1983; 1985) by employing case (B) and model (2) are the following: Table 2 shows the results of the calculations with model (2), using k(Tn-2) at x = L. The Biot number (Bif) obtained was Bimin = 17.82, with the mean value of the Biot 7.41. The typical runtime for the program is 3 s.
CR IP T
numbers of the sample Bim = 25.56, Bimax = 35.63, and the standard deviation σn-1 = The (tcalc) with or without the use of dT(1)/dt at x = 0 when using model (2) are similar or equal. Therefore, it is generally unnecessary to use dT(1)/dt at x = 0.
The (tcalc) from model (2) using dT(n)/dt at x = L is generally different from the value calculated without this derivative. Therefore, it is important to use dT(n)/dt at x = L
AN US
in model (2).
Generally, the oscillations are smaller with k((1/2)(Tn-1+Tn-2)) at x = L in model (2), than those with k(Tn-2).
The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
M
commercial calculation of (tcalc) agrees well with the experimental data (texper).
ED
4.2. Simulations using experimental data from Bazan and Mascheroni
The highlights of the prediction of (tcalc) using data from Bazan and Mascheroni (1984)
PT
by employing case (B) are as follows: Table 2 shows the results of the calculations with model (2) using k(Tn-2) at x = L.
CE
The Bif obtained were Bimin = 1.55, Bim = 1.96, Bimax = 2.35, and σn-1 = 0.33. The typical runtime is 3 s.
AC
The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or commercial calculation of (tcalc) agrees well with the experimental data (texper).
4.3. Simulations using experimental data from Cleland and Earle
The highlights of the prediction of (tcalc) using the data reported by Cleland and Earle (1977a; p. 1392) by employing case (B) are as follows:
ACCEPTED MANUSCRIPT The Bif obtained for tylose gel slabs were Bimin = 0.22, Bim = 3.56, Bimax = 20.78, and σn-1 = 5.69. The typical runtime is 4 s. The Bif obtained for mashed potato slabs were Bimin = 0.22, Bim = 1.81, Bimax = 4.21, and σn-1 = 1.59. The typical runtime is 3 s. The Bif obtained for minced lean beef slabs were Bimin = 0.27, Bim = 2.83, Bimax = 10.29, and σn-1 = 3.80. The typical runtime is 3 s. The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
CR IP T
commercial calculation of (tcalc) agrees well with the experimental data (texper).
4.4. Simulations using experimental data from Hense and Kieckbusch
Kieckbusch (1991) are the following:
AN US
The relevant aspects of the prediction of (tcalc) using the data from Hense and
The Bif obtained were Bimin = 4.05, Bim = 5.10, Bimax = 6.28, and σn-1 = 1.03. The typical runtime is 3 s.
The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
M
commercial calculation of (tcalc) agrees well with the experimental data (texper).
ED
4.5. Simulations using experimental data from Hung and Thompson
The relevant aspects of the prediction of (tcalc) using the data from Hung and Thompson
PT
(1983) are the following:
Only one experiment on tylose gel slabs, with Bif = 0.09 and a water content of
CE
0.751, required a large n and the runtime in case (B), using the data from Hung and Thompson (1983).
AC
The Bif obtained for tylose gel slabs were Bimin = 0.08, Bim = 1.15, Bimax = 3.37, and σn-1 = 0.99. The typical runtime is 5 s.
The Bif obtained for lean beef slabs were Bimin = 0.08, Bim = 1.28, Bimax = 3.30, and σn-1 = 1.09. The typical runtime is 3 s. The Bif obtained for ground beef slabs were Bimin = 0.06, Bim = 1.17, Bimax = 3.10, and σn-1 = 1.04. The typical runtime is 3 s. The Bif obtained for carp slabs were Bimin = 0.06, Bim = 1.16, Bimax = 3.15, and σn-1 = 1.08. The typical runtime is 3 s.
ACCEPTED MANUSCRIPT The Bif obtained for mashed potato slabs were Bimin = 0.06, Bim = 1.11, Bimax = 3.03, and σn-1 = 1.03. The typical runtime is 3 s. The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or commercial calculation of (tcalc) agrees well with the experimental data (texper).
4.6. Simulations using experimental data from Lescano
CR IP T
The relevant aspects of the prediction of (tcalc) using the data from Lescano (1973) are the following:
The Bif obtained were Bimin = 0.86, Bim = 3.06, Bimax = 6.18, and σn-1 = 1.59. The typical runtime is 3 s.
The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
AN US
commercial calculation of (tcalc) agrees well with the experimental data (texper).
4.7. Simulations using experimental data from Mascheroni and Calvelo
Calvelo (1982) are the following:
M
The relevant aspects of the prediction of (tcalc) using the data from Mascheroni and
ED
The Bif obtained for lean beef slabs were Bimin = 2.82, Bim = 5.62, Bimax = 8.76, and σn-1 = 1.65. The typical runtime is 3 s. The Bif obtained for lean beef slabs with heat transfer parallel to the fibers were Bimin
PT
= 1.39, Bim = 4.70, Bimax = 7.36, and σn-1 = 2.66. The typical runtime is 3 s. The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
CE
commercial calculation of (tcalc) agrees well with the experimental data (texper).
AC
4.8. Simulations using experimental data from Michelis and Calvelo
The relevant aspects of the prediction of (tcalc) using the data from Michelis and Calvelo (1983) are the following: The Bif obtained were Bimin = 2.81, Bim = 5.18, Bimax = 7.36, and σn-1 = 1.97. The typical runtime is 3 s. The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or commercial calculation of (tcalc) agrees well with the experimental data (texper).
ACCEPTED MANUSCRIPT
4.9. Simulations using experimental data from Pham and Willix
The highlights for the prediction of (tcalc) using the data from Pham and Willix (1990) are as follows: The Bif obtained were Bimin = 1.17, Bim = 7.57, Bimax = 31.93, and σn-1 = 6.59. The typical runtime is 4 s.
CR IP T
The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
commercial calculation of (tcalc) agrees well with the experimental data (texper).
5. Conclusions
AN US
A numerical model has been developed for predicting the freezing process of foods or analogues. The finite difference method (FDM) was employed to calculate the temperature distribution in a food T(x,t). The freezing times predicted (tcalc) with model (2) using the properties of tylose gel were compared with published experimental data
M
(texper).
In summary, the main conclusions of this study are as follows:
ED
It is important to use the derivative dT(n)/dt at x = L in model (2) in relation to the use of dT(1)/dt at x = 0 in the same model. The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
PT
commercial calculation of (tcalc) is adequate. Simulations using k(Tn-2) at x = L with model (2) and case (B) were performed
CE
obtaining (tcalc) and compared with 227 experimental data points (texper), exhibiting good agreement. The parameters obtained for 227 experimental data points were Emin
AC
= -6.97 %, Emean = 0.06 %, Emax = 8.08 %, and n-1 = 2.19 %.
REFERENCES
Albasiny, E.L., 1956. The solution of non-linear heat-conduction problems on the Pilot Ace. Proceedings Inst. Electric Eng. 103, Part B, issue 1, 158-162. Bazan, H.C., Mascheroni, R.H., 1984. Heat transfer with simultaneous change of phase in freezing boned mutton. Lat. Am. J. Heat Mass Transfer. 8, 55-76.
ACCEPTED MANUSCRIPT Bonacina, C., Comini, G., 1973. On the solution of nonlinear heat conduction equations by numerical methods. Int. J. Heat Mass Transfer. 16, 581-589. Carslaw, H.S., Jaeger, J.C., 1959. Conduction of Heat in Solids, Clarendon Press. Cleland, A.C., 1977. Heat transfer during freezing of foods and prediction of freezing times. Ph. D. thesis, Massey University, New Zealand. Cleland, A.C., Earle, R.L., 1977a. A comparison of analytical and numerical methods of predicting the freezing times of foods. J. Food Sci. 42(5), 1390-1395.
CR IP T
Cleland, A.C., Earle, R.L., 1977b. The third kind of boundary condition in numerical freezing calculations. Int. J. Heat Mass Transfer. 20, 1029-1034.
Cleland, A.C., Earle, R.L., 1979. A comparison of methods for predicting the freezing times of cylindrical and spherical foodstuffs. J. Food Sci. 44, 958-963, 970.
Cleland, A.C., Earle, R.L., 1984. Assessment of freezing time prediction methods. J.
AN US
Food Sci. 49, 1034-1042.
Cleland, D.J., Cleland, A.C., Jones, R.S., 1994. Collection of accurate experimental data for testing the performance of simple method for food freezing time prediction. J. Food Proc. Eng. 17, 93-119.
M
Crank, J., Nicolson, P., 1947. A practical method for the numerical integration of solutions of partial differential equations of heat conduction type, Proceedings Cambridge Philosophical Society, 43, 50–67.
ED
Creed, P.G., James, S.J., 1983. The freezing times of liver in a vertical plate freezer. Progress in Refrigeration Science and Technology, Proceedings of the 16th
PT
International Congress of Refrigeration, Tome III, Paris, France, 4, 145-151. Creed, P.G., James, S.J., 1985. Heat transfer during the freezing of liver in a plate
CE
freezer. J. Food Sci. 50, 285-288, 294. Cutlip, M.B., Shacham, M., 1999. Problem Solving in Chemical Engineering with
AC
Numerical Methods, Prentice Hall. Dehghan, M., 2006. A computational study of the one-dimensional parabolic equation subject to nonclassical boundary specifications. Num. Meth. Partial Dif. Eq. 22, 220257.
Dima, J.B., Santos, M.V., Baron, P.J., Califano, A., Zaritzky, N.E., 2014. Experimental study and numerical modeling of the freezing process of marine products. Food Bioproducts Processing. 92, 54-66. Ferreira, S.R., Rojas, L.O.A., Souza, D.F.S., Oliveira, J.A., 2016. Freezing time of an infinite cylinder and sphere using the method of lines. Int. J. Refrigeration. 68, 37-49.
ACCEPTED MANUSCRIPT Ferreira, S.R., 2017. Freezing time of a slab using the method of lines. Int. J. Refrigeration. 75, 77-94. Hense, H., Kieckbusch, T.G., 1991. Freezing of shark: I Experimental results. In Proceedings of the IV Latin American Congress of Heat and Mass Transfer. 112 Hung, Y.C., Thompson, D.R., 1983. Freezing time prediction for slab shape foodstuff by an improved analytical method. J. Food Sci. 48, 555-560.
Corporate Headquarters, Houston, Texas.
CR IP T
Imsl, 2006. Fortran numerical library user‟s guide, version 6.0, Visual Numerics
Lees, M., 1966. A linear three level difference scheme for quasi-parabolic equations. Maths. Comput. 20, 516–522.
Lescano, C.E., 1973. Predicting freezing curves in codfish fillets using the ideal binary solution assumption. MS thesis, Michigan State University.
AN US
Luikov, A.V., 1968. Analytical Heat Diffusion Theory, Academic Press.
Mannapperuma, J.D., Singh, R.P., 1988. Prediction o freezing and thawing times of foods using a numerical method based on enthalpy formulation. J. Food Sci. 53, 626630.
M
Mannapperuma, J.D., Singh, R.P., 1989. A computer-aided method for prediction of properties and freezing/thawing times of foods. J. Food Eng. 9, 275-304. Mascheroni, R.H., Calvelo, A., 1982. A simplified model for freezing time calculations
ED
in foods. J. Food Sci. 47, 1201-1207. Michelis, A., Calvelo, A., 1983. Freezing time predictions for brick and cylindrical-
PT
shaped foods.J. Food Sci. 48, 909-913, 934. Northrop, P.W.C., Ramachandran, P.A., Schiesser, W.E., Subramanian, V.R., 2013. A
CE
robust false transient method of lines for elliptic partial differential equations. Chem. Eng. Sci. 90, 32-39.
AC
Pham, Q.T., 1985. A fast, unconditionally stable finite-difference scheme for heat conduction with phase change. Int. J. Heat Mass Transfer. 28(11), 2079-2084.
Pham, Q.T., 2008. Modelling of freezing processes, in: Evans, J.A. (Ed.), Frozen Food Science and Technology. Blackwell, Oxford, U.K., pp. 51-80.
Pham, Q.T., 2014. Food Freezing and Thawing Calculations, Springer, New York. Pham, Q.T., Willix, J., 1990. Effect of Biot number and freezing rate on accuracy of some food freezing time prediction methods. J. Food Sci. 55(5), 1429-1434. Sadiku, M.N.O., Obiozor, C.N., 2000. A simple introduction of the method of lines. Int. J. Electrical Eng. Educ. 37(3), 282-296.
ACCEPTED MANUSCRIPT Salvadori, V.O., 1994. Heat transfer during freezing, storage and thawing of foods. Ph. D. thesis, University National of La Plata, Argentina. (In Spanish.) Schiesser, W.E., 1991. The Numerical Method of Lines: Integration of Partial Differential Equations, Academic Press, San Diego. Schiesser, W.E., Silebi, C.A., 2009. Computational Transport Phenomena: Numerical Methods for the Solution of Transport Problems, Cambridge University Press. Schwartzberg, H.G., Singh, R.P., Sarkar, A., 2007. Freezing and thawing of foods –
CR IP T
computation methods and thermal properties correlation, in: Yanniotis, S., Sunden, B. (Eds.), Heat Transfer in Food Processing: Recent Developments and Applications, WIT Press, pp. 61-99.
Singh, R.P., Heldman, D.R., 2009. Introduction to Food Engineering, Elsevier.
Succar, J., Hayakawa, K., 1983. Empirical formulae for predicting thermal physical
AN US
properties of food at freezing or defrosting temperatures. Lebensm.-Wiss.u-Technol. 16(6), 326-331.
Wilson, H.A., Singh, R.P., 1987. Numerical simulation of individual quick freezing
AC
CE
PT
ED
M
spherical foods. Int. J. Refrigeration. 10, 149-155.
ACCEPTED MANUSCRIPT
i =1 T(1,t)
2 T(2,t)
3 T(3,t)
(i-1) T(i-1,t)
i T(i,t)
(i+1) T(i+1,t)
CR IP T
Fig. 1 – Grid of the half-thickness of a slab.
Table 1 – Computer program and subroutine for model (2) using tylose gel properties by a selected experiment. Program SlabDerivativeCaseB
M
PT
ED
Write(10,10) tend, (y(i), i = 1, n) t = tend tend = t + step End Do 10 Format ((2x, f8.2)) !Format to print results. End program Subroutine Fcn(n, t, y, yprime) Implicit none Integer i, j, n Real*8 t, Y(n), Yprime(n) Real*8 k, Ya, Ym, h, L, dx, Rho, Cp Ya = -21 h = 361 L = 0.0380 dx = (L - 0.d0)/dble(n - 1)
AN US
Include 'link_fnl_shared.h' Use IVMRK_INT Implicit none Integer i, j, ido, n, np Parameter (n = 7, np = 5000) real*8 t, tend, Y(n), Yprime(n), step External Fcn ! Set initial values. t = 0.d0 Y(1:n) = 10.5 ! Defining final time and integration step. tend = 10000 step = (tend -t)/dble(np) tend = step ! Creating an output data file. Open (10, file = "Results.txt", status = "unknown") Write(10,10) t, (y(i), i = 1, n) Ido = 1 ! Normally, the initial call is made with IDO = 1. Do j = 1, np Call DIvmrk(ido, n, fcn, t, tend, y, yprime)
CE
!Ym = T(average) selected for calculation in x = L and i = n. !Ym = Y(n-1) Ym = Y(n-2) !Ym = (1/2)*(Y(n-2) + Y(n-1)) !Ym = (1/3)*(Y(n-2) + Y(n-1) + Y(n))
AC
If (Y(1).GT.-0.6) Then k = 0.49 Rho = 1006.0 Cp = 3688.0 Else If (Y(1).LE.-0.6) Then k = 1.573 - 0.00773*Y(1) + 0.653/Y(1) Rho = 937.908 - 0.154*Y(1) - 40.80/Y(1) Cp = 1418 + 122173/((-Y(1))**1.696) End If ! Evaluate derivative of temperature in x = 0 and i = 1: !Y(1) = Y(2) Yprime(1) = k*(Y(2)- 2*Y(2) + Y(3))/(Rho*Cp*(dx**2)) Do i = 2, n - 1 If (Y(i).GT.-0.6) Then k = 0.49 Rho = 1006.0 Cp = 3688.0 Else If (Y(i).LE.-0.6) Then k = 1.573 - 0.00773*Y(i) + 0.653/Y(i) Rho = 937.908 - 0.154*Y(i) - 40.80/Y(i) Cp = 1418 + 122173/((-Y(i))**1.696) End If
i=n T(n,t)
ACCEPTED MANUSCRIPT
! Evaluate derivatives of temperature dT/dt: Yprime(i) = k*(Y(i-1) - 2*Y(i) + Y(i+1))/(Rho*Cp*(dx**2)) End Do If (Ym.GT.-0.6) Then k = 0.49 Rho = 1006.0 Cp = 3688.0 Else If (Ym.LE.-0.6) Then k = 1.573 - 0.00773*Ym + 0.653/Ym Rho = 937.908 - 0.154*Ym - 40.80/Ym Cp = 1418 + 122173/((-Ym)**1.696) End If
CR IP T
! Evaluate derivative of temperature in x = L and i = n: !k*(Y(n)-Y(n-1)) = dx*h*(Ya-Y(n)) Yprime(n)= (k*(Y(n-2)-Y(n-1)) + dx*h*(Ya-Y(n)))/(Rho*Cp*(dx**2.0)) End subroutine
for slabs using case (B) and model (2).
AN US
Table 2 – Summary of comparisons between experimental and predicted freezing times
Material and model
Evaluation of (k) in x = L
Tc (°C)
Runs
Xmin (%)
Xmean (%)
Xmax (%)
n-1 (%)
Bazan and Mascheroni (1984)
Boned mutton
k(Tn-2)
-18
4
-1.25
-0.21
0.22
0.70
Cleland and Earle (1977a)
Tylose gel
k(Tn-2)
-10
43
-3.44
0.89
4.98
2.08
Cleland and Earle (1977a)
Mashed potato
k(Tn-2)
-10
6
-2.05
0.00
3.21
1.78
Cleland and Earle (1977a)
Minced lean beef
k(Tn-2)
-10
6
-1.98
0.06
2.31
1.68
Creed and James (1983; 1985)
Pig liver
k(Tn-2)
-7
44
-3.87
0.00
5.14
1.91
Pig liver
k((1/2)(Tn-2+Tn-1))
-7
44
-2.94
-0.20
1.69
0.76
Shark meat
k(Tn-2)
-10
6
-3.44
-0.83
1.33
1.71
ED
PT
CE
Hense and Kieckbusch (1991)
M
Source of data
Tylose gel
k(Tn-2)
-18
23
-6.97
-0.79
2.24
2.00
Hung and Thompson (1983)
Lean beef
k(Tn-2)
-18
9
-2.09
0.34
2.23
1.24
Hung and Thompson (1983)
Ground beef
k(Tn-2)
-18
9
-4.35
-0.89
5.24
2.70
Hung and Thompson (1983)
Carp
k(Tn-2)
-18
9
-5.37
-0.59
2.01
2.29
Hung and Thompson (1983)
Mashed potato
k(Tn-2)
-18
9
-4.20
-0.37
2.22
2.18
Lescano (1973)
Codfish fillet
k(Tn-2)
-18
8
-5.65
0.42
2.34
2.72
Mascheroni and Calvelo (1982)
Lean beef
k(Tn-2)
-18
10
-3.11
0.22
3.34
2.38
AC
Hung and Thompson (1983)
ACCEPTED MANUSCRIPT
Lean beef (heat transfer parallel to the fibers)
k(Tn-2)
-18
4
-3.82
0.50
8.08
5.48
Michelis and Calvelo (1983)
Lean beef
k(Tn-2)
-18
5
-3.10
-0.13
3.34
2.98
Pham and Willix (1990)
Tylose gel
k(Tn-2)
-18
32
-3.81
0.20
5.53
2.28
AC
CE
PT
ED
M
AN US
CR IP T
Mascheroni and Calvelo (1982)
Freezing times using time derivative of temperature on surface of foods S.R. Ferreira , L.O.A. Rojas PII: DOI: Reference:
S0140-7007(18)30442-0 https://doi.org/10.1016/j.ijrefrig.2018.11.009 JIJR 4168
To appear in:
International Journal of Refrigeration
Received date: Revised date: Accepted date:
2 June 2018 26 October 2018 7 November 2018
Please cite this article as: S.R. Ferreira , L.O.A. Rojas , Freezing times using time derivative of temperature on surface of foods, International Journal of Refrigeration (2018), doi: https://doi.org/10.1016/j.ijrefrig.2018.11.009
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Highlights A numerical model was developed to predict the food freezing process. (FDM) was employed to calculate the temperature distribution in a food T(x,t). This paper tests the use of dT(i=n)/dt in x=L, and dT(i=1)/dt in x=0. Simulations obtaining (tcalc) were compared with 227 experimental data (texper).
AC
CE
PT
ED
M
AN US
CR IP T
Good agreement was obtained between predicted and experimental values.
ACCEPTED MANUSCRIPT Freezing times using time derivative of temperature on surface of foods S.R. Ferreira a*, L.O.A. Rojas b a
Department of Chemical Engineering, UFRN - University Federal of the Rio Grande
do Norte, Campus Universitario - 59078-970 - Natal - RN - Brazil. Tel. 55-8432153753 b
Department of Chemical Engineering, UFPb – Federal University of the Paraiba -
*
[email protected]
Abstract
CR IP T
Joao Pessoa - Pb - Brazil.
A numerical model is developed to predict the freezing process of foods or analogues.
AN US
The finite difference method (FDM) is employed to calculate the temperature
distribution in a food T(x,t). The aim of this research study is to test the use of the time derivative of temperature on the surface, dT(i=n)/dt in x=L, and at the center of food, dT(i=1)/dt in x = 0. This methodology is employed to predict the freezing time (tcalc),
M
using the properties of tylose gel (k, ρ, and Cp), which is a food analog. Simulations were performed using the model, obtaining (tcalc), and compared with 227 experimental data points (texper). In summary, the main conclusions of this study are as follows. It is
ED
important to use the derivative dT(n)/dt at x = L in model (2) in relation to the use of dT(1)/dt at x = 0 in the same model. The evaluation of thermal conductivity using
PT
temperature T = Tn-2, k(Tn-2) at x = L with the proposed model and properties of tylose gel to calculate (tcalc) is adequate. Simulations using k(Tn-2) at x = L with the proposed
CE
model were performed, obtaining (tcalc), and compared with 227 experimental data points (texper), exhibiting good agreement. The following parameters were obtained for
AC
227 experimental data points: minimum error Emin = -6.97 %, overall mean error Emean = 0.06 %, maximum error Emax = 8.08 %, and standard deviation n-1 = 2.19 %. Keywords: time derivative of temperature; convective boundary condition; freezing time; finite difference method.
Nomenclature
Bi
Biot number [--]
ACCEPTED MANUSCRIPT Bif
Biot number in the fully frozen state [--]
Bim
mean value of Biot numbers of the sample in the fully frozen state [--]
Bimax maximum value of Biot number in the fully frozen state [--] minimum value of Biot number in the fully frozen state [--]
Cp
effective specific heat [J kg-1 °C-1]
Cp0
unfrozen effective specific heat [J kg-1 °C-1]
Cp1
freezing zone effective specific heat [J kg-1 °C-1]
dx
step size for spatial variable [m]
E
error in calculation or experimental percentage error [%]
h
convective heat transfer coefficient [W m-2 °C-1]
H
specific enthalpy [J kg-1]
i
ith node [--]
k
thermal conductivity [W m-1 °C-1]
kf
thermal conductivity of food in fully frozen state [W m-1 °C-1]
k(i)
thermal conductivity calculated at temperature of ith node T(i) [W m-1 °C-1]
k0
unfrozen thermal conductivity [W m-1 °C-1]
k1
freezing zone thermal conductivity [W m-1 °C-1]
L
slab half-thickness [m]
n
number of nodal points or on surface [--]
np
number of divisions in time [--]
R
radius of infinite cylinder or sphere, or object half-thickness (L) [m]
t
time [s]
tcalc
calculated freezing time [s]
texper
experimental freezing time [s]
AN US
M
ED
PT
integration time [s] temperature [°C]
AC
T
CE
tend
CR IP T
Bimin
Ta
ambient or cooling medium temperature [°C]
Tc
center temperature or desired final (T) at center of food [°C]
Tf
initial freezing temperature [°C]
Tm
average temperature at surface in i = n for k(Tm) calculation [°C]
Tn
temperature at surface with convection [°C]
T0
initial temperature [°C]
x
distance measured along the x-axis [m]
Xmax
maximum deviation [%]
ACCEPTED MANUSCRIPT Xmean mean deviation [%] Xmin
minimum deviation [%]
Y
temperature in computational program [°C]
Ya
cooling medium temperature in computational program [°C]
Ym
average temperature at surface in i = n for k(Tm) calculation [°C]
Greek symbols
distance increment in (x) dimension [m]
ρ
density [kg m-3]
ρ0
unfrozen density [kg m-3]
ρ1
freezing zone density [kg m-3]
σn-1
standard deviation [%] or [--]
M
Subscripts
AN US
∆x
ambient
c
center
calc
calculated value
end
end time of calculation in computational program
ED
a
PT
exper experimental value freezing
i
ith
m
mean or average maximum
AC
max
CE
f
min
minimum
n
number of nodes in (x) dimension
0
initial or unfrozen
1
in freezing zone
1. Introduction
CR IP T
Yprime time derivative T/t in computational program [°C s-1]
ACCEPTED MANUSCRIPT Freezing times (tcalc) are estimated to predict the required residence times for foods in freezers. As foods freeze, a frozen layer forms next to the outer surface and moves progressively inward. Ice grows, and heat is removed, mainly at, or near the inwardmoving freezing interface (Schwartzberg et al., 2007, p. 70). Commercial finite element method (FEM) software often provide graphical user interfaces, are user-friendly, can be used without understanding FEM principles, and have been used to treat complex and noncomplex freezing and thawing problems
CR IP T
(Schwartzberg et al., 2007, p. 89).
Several numerical methods have been developed to solve partial differential
equations (Northrop et al., 2013). For parabolic equations such as the heat equation,
there exist several numerical methods to find a solution (Dehghan, 2006; Northrop et
al., 2013). The method of lines (MOL) is one such efficient routine in which the space is
AN US
discretized using any of a number of techniques such as the finite difference method (FDM), finite element method (FEM), finite volume method (FVM), or collocation method (CM) (Cutlip and Shacham, 1999; Dehghan, 2006; Northrop et al., 2013;
Sadiku and Obiozor, 2000; Schiesser, 1991; Schiesser and Silebi, 2009). Ferreira et al.
M
(2016) and Ferreira (2017) used the MOL to predict (tcalc), based on the effective specific heat method.
Generally, commercial software based on the FDM are unable to run and/or converge
ED
to accurate temperature predictions when attempting to simulate a phase change process, owing to abrupt changes in the thermophysical properties, specifically when
PT
using the apparent specific heat method (Dima et al., 2014). The surface of an object is the most difficult to model in the freezing process when
CE
there is a convective boundary condition (Albasiny, 1956; Bonacina and Comini, 1973, p. 586-587; Cleland, 1977, p. 18-19; Cleland and Earle, 1977b, p. 1030-1032). For this
AC
reason, the focus of this approach is to determine the influence of the time derivatives of temperature, dT(i=n)/dt and dT(i=0)/dt. This approach allows the calculation of the heat transfer convective coefficient (h) and an accurate prediction of the (tcalc). Based on what has been cited, the objectives of this research are as follows: To develop a mathematical model, called model (2), for simulating the freezing process in a slab, and calculate the (tcalc). To evaluate the influence of the use of dT(i=n)/dt and dT(i=0)/dt. To compare the predictions of the proposed model (tcalc) to published experimental data (texper).
ACCEPTED MANUSCRIPT Experimental freezing times were selected from published data (Bazan and Mascheroni, 1984; Cleland and Earle, 1977a; Creed and James, 1983, 1985; Hense and Kieckbusch, 1991, 1985; Hung and Thompson, 1983; Lescano, 1973; Mascheroni and Calvelo, 1982; Michelis and Calvelo, 1983; Pham and Willix, 1990).
2. Literature review
CR IP T
In this section, some published numerical methods that are used to predict (tcalc), and several challenges of the phase change process modeling are briefly reviewed.
2.1. Some numerical methods for freezing-time prediction
AN US
A survey of the literature since 1984 shows that this systematic approach for assessing prediction methods has not been widely used (Cleland et al., 1994). This survey also shows that no method (including numerical methods) predicts all experimental data well, indicating that some data that are widely used for testing are of lower quality than
M
that of others. Furthermore, it may be possible to ensure that the data collected in the future is more accurate (Cleland et al., 1994). The choice of technique depends on the objectives of the modelers and the technical
ED
means at their disposal (Pham, 2008, 2014). Numerical methods can, in principle, provide exact or near-exact predictions for almost any situation, although in practice,
PT
their accuracy is limited by an inadequate knowledge of the problem parameters, for example, the product properties, geometry, and flow characteristics (Pham, 2008, 2014).
CE
Schwartzberg et al. (2007) discussed the apparent specific heat formulation and the enthalpy step method for freezing and thawing calculations, including the numerical
AC
solution of the partial differential equations (PDEs) describing the heat transfer during these processes. One such method, based on the use of specific heat (Cp), calculates the T(x,t) using
explicit methods, such as Euler‟s forward difference method. The following information provided by Schwartzberg et al. (2007) is relevant for this method of calculating the (tcalc): The (tcalc) and T(x,t) for food freezing, obtained using these methods, are presented.
ACCEPTED MANUSCRIPT The computation becomes unstable if the time step (t) is too large. Stability is provided by using the smallest (t) obtained from the equations presented in the literature (Mannapperuma and Singh, 1988, 1989; Schwartzberg et al., 2007). When the equations for obtaining T(x,t) were solved for typical freezing conditions by machine computation, 700–2,000 time steps were used for slabs, 370–1,000 steps for infinite cylinders, and 130–1,500 steps for spheres. The number of steps needed is inversely proportional to (L2) and proportional to
CR IP T
n2(n+1), but also depends on Bif = 2hL/kf, (T0), (Ta), and (Tc); (n) is the number of
nodal points for FDM calculations, (Bif) is Biot number in the fully frozen state, (h) the convective heat transfer coefficient, (L) the object‟s half thickness or radius (R), (kf) the thermal conductivity of the food in the fully frozen state, (T0) the initial
AN US
temperature, (Ta) the temperature of the cooling, and (Tc) the desired final (T) at the center of the food. The typical Biot number is defined as Bi = hL/k (Luikov, 1968; Singh and Heldman, 2009); which is different from that employed by Schwartzberg et al. (2007) and Pham and Willix (1990).
One can use a much larger (t) with both the Lee and Crank–Nicolson methods than
M
that for the explicit method, but if (t) is too large, „jumping‟ and (T) oscillations
ED
occur as (T) passes through (Tf), where (Tf) is the initial freezing point.
2.2. Challenges of the phase change process modeling
PT
The problem of predicting the time for a material to freeze under a given set of conditions is classed physically as one of heat conduction with change of phase
CE
(Cleland, p. 18, 1977).
The change of phase can be taken into account in two different ways. It can be
AC
assumed that all the latent heat is released at a freezing temperature at which a step change in thermophysical properties occurs (Carslaw and Jaeger, p. 282, 1959; Cleland, p. 18, 1977). Alternatively, the release of latent heat over a range of temperatures can be taken into account by changing specific heat Cp(T) and thermal conductivity k(T). The first has been used extensively, even though variations in thermal properties mean that a freezing foodstuff will fit the solution only approximately. The second method introduces a problem in defining the completion of freezing, which is deemed complete when the phase change front reaches the thermodynamic centre of the material in the
ACCEPTED MANUSCRIPT first method, and the presence of a mushy zone precludes the use of such a definition (Cleland, p. 19, 1977). There are several challenges for modelers in the food freezing process, such as trying to understand the phase change, mathematical modeling of the process and trying to solve the numerical problem to calculate the freezing time. The freezing of foodstuffs frequently involves boundary conditions of the third kind, which have proved to be difficult to handle and subject to error, mainly because of the
CR IP T
requirement to calculate surface temperatures implicitly (Cleland and Earle, p. 1029, 1977b).
The focus of this section is the boundary condition on the surface of a slab, because it is the most difficult point to be modeled, that is, to analyze the use of derivative dT(i=n)/dt. The use of derivative dT(i=1)/dt in x = 0 is also tested.
AN US
In design situations the surface temperature/time profile is rarely known. However
the bulk average ambient temperature (Ta) is often known, for example, air temperature in a blast freezer. There are also well known means of arriving at a surface heat-transfer coefficient (h), which considers all resistance for heat leaving the surface of the freezing
M
material and entering the cooling medium (Cleland and Earle, p. 1030, 1977b). To facilitate analysis, a slab is considered. The half-thickness of the slab (L) in the range of 0 ≤ x ≤ +L is divided into (n-1) sections with (n) node points, numbered from
ED
(1) to (n); that is, setting i = 1 in the center of the slab when x = 0, and i = n on the slab surface when x = L. The space grid is shown in Fig. 1.
PT
In these cases, the surface requires special consideration. To solve some of the difficulties and improve freezing time accuracy, a number of researchers proposed
CE
changes in the mathematical modeling on the surface of a slab. Some of these are described below:
AC
Bonacina and Comini (1973, p. 586) suggest that when a fictitious point i = (n+1) is used, then T(n+1) = Ta, and pseudo-thermal conductivity is defined, k(n+1/2) = hx, so that a three level scheme can be used. The authors found the “Lees” scheme to be superior to others, due to smaller truncation errors, and better representation of thermal properties (Cleland and Earle, 1977b, p. 1030). It is important to underscore, however, that the three-time level method remains unconditionally stable, even with boundary conditions of the second and third kind (Bonacina and Comini, 1973, p. 587).
ACCEPTED MANUSCRIPT In another approach, the finite difference boundary condition was derived by Cleland and Earle (1977b, p. 1032) from a heat balance over the surface space increment which extends (0.5x) from the surface. While the internal nodes have the specific heat (Cp) of a full (x) associated with them, the surface node has the (Cp) of only a half space increment. The freezing times calculated from this approach were compared to their experimental data (Cleland, 1977, p. 103).
CR IP T
Albasiny (1956) used the enthalpy method and a fictitious point i = (n+1) with T(n+1) = Ta, obtaining a time derivative of enthalpy on the surface of the slab,
dH(i=n)/dt in x=L. The set of resulting equations was solved using the Runge-Kutta method, obtaining enthalpy (H) and the corresponding temperature (T) versus time (t). According to Cleland (1977, p.97) the solution obtained by Albasiny (1956) is
due to the large truncation error.
3. Materials and methods
AN US
not as accurate as the scheme of Lees (1966) or the Crank-Nicolson (1947) method
M
The method is described in five sections. The freezing process model is described in Section 3.1, based on the heat conduction equation for a slab with the denominated
ED
model (2). This model is discretized using an FDM, as shown in Section 3.2, and solved using a computational program by the MOL, as shown in Section 3.3. The selected experimental (tcalc) dataset from the literature is presented in Section
PT
3.4. Finally, the selection of the physical properties that are employed as parameters in the numerical program when implementing model (2) is discussed in Section 3.5. The
CE
results obtained are discussed in the Results and discussion section.
AC
3.1. Unsteady-state freezing for symmetric heat conduction in slabs
For illustrative purposes, the freezing process in slab-shaped foods is analyzed. Analogous equations can be written for spheres and infinite cylinders. A slab located within the range of -L ≤ x ≤ +L is considered. When the heat transfer is symmetric, it is convenient to use the slab center as the origin for the coordinates. Because of symmetry, only the region between x = 0 and x = L is treated. The half-
ACCEPTED MANUSCRIPT thickness of the slab (L) considered is divided into (n-1) sections with (n) node points, indexed from (1) to (n) from the center to the slab border. A heat balance over a small section of a slab of thermal conductivity k(T), assuming that the material within the ith element has a specific heat (Cp) and density (ρ), can be represented by (Cleland, 1977; Cleland and Earle, 1984)
()
()
( ( )
⁄
)
⁄
At the limit when ∆x 0, Eq. (1) gives
()
()
( ( )
)
⁄
)
⁄
(1)
(2)
AN US
()
( ( )
CR IP T
()
The first term of Eq. (1) or (2), that is, the accumulation of heat in the slab, can be calculated at a temperature (Ti) with Cp(Ti) and (Ti). The last term, representing the heat conduction in the slab, can be evaluated at positions i+1/2 and i−1/2. There are
M
various options for choosing the (T), for which the properties (k), (), and (Cp) are calculated. For example, the second term of Eq. (2) can be calculated at the “centered”
ED
position (i), k(i) = k(Ti).
For a slab with thickness (2L) that is cooled from both sides at a cooling temperature
PT
(Ta), the initial and boundary conditions are
CE
(3)
AC
( )
( )
(
(4)
)
(5)
Model (2) is discussed as follows. The importance of dT(i=1)/dt and dT(i=n)/dt for the freezing-process modeling is tested. The procedure used is as follows: The slab in the range 0
x
L of Eq. (2) is discretized; this is valid from i = 1 to i =
n. However, it is usually discretized only from i = 2 to i = n-1, using the boundary
ACCEPTED MANUSCRIPT conditions at the extremes, i = 1 and i = n, to complete the energy balance for the slab at 0
x
L.
In this paper, Eq. (2) is also discretized for i = 1 and x = 0. The boundary condition of Eq. (4) is discretized and substituted into Eq. (2). It results in a discretized form for i = 1, with the first member equal to dT(i=1)/dt. Additionally, Eq. (2) is discretized for i = n and x = L. The boundary condition of Eq.
with the first member equal to dT(i=n)/dt.
CR IP T
(5) is discretized and substituted into Eq. (2). It results in a discretized form for i = n, With this procedure, all the interior and extreme points of the discretized equations will have a "local equation" as a function of dT(1)/dt, dT(2)/dt, …, dT(n-1)/dt, and dT(n)/dt.
AN US
The resulting system of equations is solved using the initial and boundary conditions. Solving the system yields the temperatures T(i) versus time (t), that is, T(1,t), T(2,t), …, T(n-1,t), and T(n,t).
M
3.2. Discretization of model (2) for slabs
After realizing the discretizations for the space, the resulting system of ordinary
ED
differential equations (ODE) is solved by the MOL. The half-thickness of the slab (L) in the range of 0 ≤ x ≤ +L is divided into (n-1)
CE
given by
PT
sections with (n) node points, numbered from (1) to (n); therefore, the space grid is
(6)
)
AC
(
Setting i = 1 in the center of the slab when x = 0, and i = n on the slab surface when x
= L, Eq. (3) can be written as ()
( )
(7)
Eq. (2) can be discretized using a central second-order approximation around point (i). The properties ρ(i), Cp(i), and k = k(i) are evaluated at T(i) for the interior points 2 ≤ i ≤ n-1, resulting in
ACCEPTED MANUSCRIPT
( )(
()
) ()
( )(
(8)
)
Salvadori (1994) presents a similar space discretization for a slab, although they performed a dual time and space discretization. However, the temperature (T) at which (k) is calculated varies from author to author; for example, Cleland (1977, p. 108) and
CR IP T
Cleland and Earle (1979, p. 962) used k(i) = (1/2)(k(i-1/2)+k(i+1/2)), while Wilson and Singh (1987, p. 151) used k(i-1/2) = k((1/2)(T(i)+T(i-1))) and k(i+1/2) = k((1/2)(T(i)+T(i+1))) etc.
(
)
(
)( ( ) ( )(
) )
AN US
Setting i = 1 in the center of the slab when x = 0, Eq. (2) can be discretized as
(9)
The properties ρ(1), Cp(1), and k = k(1) were evaluated at temperature T(1) using Eq. (9). The subsequent evaluation of the properties ρ(1), Cp(1), and k = k(2) could also be
M
tested, considering that (k) is evaluated at a mean temperature T(2) etc. Various discretization schemes can be tested for the heat flow qx = -kdT/dx for x = 0
)
(10)
PT
( )( ( )
ED
using Eq. (4) (Cutlip and Shacham, 1999, p. 396). Discretizing Eq. (4) at i = 1,
CE
Rearranging Eq. (10),
AC
(11)
Substituting Eq. (11) in Eq. (9),
(
)
( )( ( )
) ( )(
)
( )( ( )
)⁄( ) ( )( )
(12)
ACCEPTED MANUSCRIPT The energy balance is complete at x = 0 when substituting Eq. (11) in Eq. (9), which results in Eq. (12), that is, T(i=1)/t. This is due to the symmetry of the temperature profile T(x,t) at i = 1, resulting in qx = 0 at i = 1, that is, T(i=1)/t at x = 0 has less influence on the calculations of T(x,t) when compared with that of T(i=n)/t on the surface of the slab at x = L. Setting i = n on the slab surface when x = L, Eq. (2) can be backward discretized as (
)
) (
(
) (
)
( )
)
CR IP T
(
(13)
The properties (, Cp, and k) can be evaluated at a selected (T) using Eq. (13). For example, (n), Cp(n), and k(n); (n), Cp(n), and k(n-1); (n-1), Cp(n-1), and k(n-2);
AN US
(n-2), Cp(n-2), and k(n-2) etc.
Several discretizations can be tested for the heat flow dT/dx at x = L using Eq. (5) (Cutlip and Shacham, 1999, p. 396). Discretizing the boundary condition of Eq. (5) at i
(
)
( )(
)
)
(14)
ED
(
M
= n at x = L gives
)
(
)( )(
)
(15)
CE
(
PT
Rearranging Eq. (14),
AC
Substituting Eq. (15) in Eq. (13), (
)
(
)
)( )(
( (
)
)
(16)
The energy balance is complete at x = L when substituting Eq. (15) in Eq. (13), which results in Eq. (16). Generally, T(i=n)/t at x = L is more important in calculations with a convective boundary condition in comparison with T(i=1)/t at x = 0, as already mentioned. To evaluate (k) and calculate () and (Cp) to use in Eq. (16), Tm = Tn-1, Tm = Tn-2, Tm
ACCEPTED MANUSCRIPT = (1/2)(Tn-2 +Tn-1) etc. are tested, where (Tm) is the average temperature in the evaluation of k(Tm) at i = n. Eqs. (6), (7), (8), (12), and (16) can be numerically integrated to any time (t) with one adequate solver. Solving the system yields temperature T(i) versus time (t), with i=1 to i=n. The number of discretization points (n) for each experiment is found by trial and error, until the lowest percentage error E(%) is obtained between the experimental
CR IP T
(texper) and calculated values (tcalc).
3.3. Development of a computational program for slabs based on model (2)
A computational program and auxiliary subroutine were implemented in Fortran 90 to integrate the resulting system of equations, as shown in Table 1, for model (2).
AN US
The program was implemented in Intel Visual Fortran, by employing the IMSL -
Fortran Numerical Library (IMSL, 2006), using routine IVMRK. All the numerical runs were tested on a Compaq Presario CQ40 Notebook PC equipped with an Intel Core 2 Duo T6500 processor with a speed of 2.1 GHz, 4 GB of RAM, and a 64-bit operating
M
system.
ED
3.4. Selected experimental freezing time (texper) dataset from the literature Certain freezing time experimental data were selected for simulation from published
PT
data (Bazan and Mascheroni, 1984; Cleland and Earle, 1977a; Creed and James, 1983, 1985; Hense and Kieckbusch, 1991; Hung and Thompson, 1983; Lescano, 1973;
CE
Mascheroni and Calvelo, 1982; Michelis and Calvelo, 1983; Pham and Willix, 1990). Several of these experimental data are cited and analyzed in the work of Salvadori
AC
(1994).
Each freezing time value that was predicted with the proposed model (tcalc) was
compared with published experimental data (texper). The results of the calculations using model (2) are discussed in the Results and discussion section.
3.5. Process of selection of physical properties used in calculations
ACCEPTED MANUSCRIPT The set of physical properties of tylose gel (k, ρ, and Cp) employed in the freezing simulation using model (2) is called case (B) (Succar and Hayakawa, 1983; Ferreira, 2017).
4. Results and discussion
CR IP T
4.1. Simulations using experimental data from Creed and James
Table 1 illustrates the computational program and subroutine for the calculation of the slab temperature T(x,t) and freezing time (tcalc) by employing model (2).
In the 44 freezing experiments on pig liver slabs conducted by Creed and James
temperature of Tc = -7 °C is reached.
AN US
(1983; 1985), the freezing time (tcalc) is the time from the onset of cooling, until a center Table 1 presents the input parameters that are necessary for calculation with the first experiment from Creed and James (1983; 1985). The parameters for the pig liver slab are n = 7, np = 5000, tend = 10000 s, L = 0.0380 m, Y(1:n) = 10.5 °C, Ya = -21 °C, and h
M
= 361 W m-2 °C-1. The time needed to reach a value of Tc = -7 °C at the center point of the pig liver slab was calculated, resulting in tcalc = 9454.0 s = 2.63 h < tend = 10000 s. The experimental time texper = 2.57 h, resulting in a percent relative error E = 2.18 %.
ED
The run time for the program is 2 s.
The criterion employed by Pham and Willix (1990, p. 1429) was used, for which the
PT
Biot number in the frozen state Bif = h(2L)/kf is considered large if Bif ≥ 10; kf is the frozen thermal conductivity. The definition Bif = h(2L)/kf is used by numerous
CE
researchers in their studies on food freezing processes (Cleland, 1977; Pham and Willix, 1990; Schwartzberg et al., 2007). Pham and Willix (1990) used a kf = 1.54 W m-1 °C-1
AC
for tylose gel. Using this value, a Bif = 17.82 ≥ 10 is obtained for the original Experiment 1 of Creed and James (1983; 1985). The values for all the other experiments are Bif ≥ 10 (Creed and James, 1983; 1985). Table 2 presents a summary of the comparisons of the experimental and predicted
freezing times for the selected case (B), using model (2). These are some of the tests performed using model (2): Analysis of the importance of dT(1)/dt at x = 0 and dT(n)/dt at x = L in the computer program of model (2).
ACCEPTED MANUSCRIPT Comparison of the use of k(Tn-2) for the evaluation of k at x = L in relation to the use of k((1/2)(Tn-2+Tn-1)) in model (2). The relevant aspects of the prediction of (tcalc) using the data from Creed and James (1983; 1985) by employing case (B) and model (2) are the following: Table 2 shows the results of the calculations with model (2), using k(Tn-2) at x = L. The Biot number (Bif) obtained was Bimin = 17.82, with the mean value of the Biot 7.41. The typical runtime for the program is 3 s.
CR IP T
numbers of the sample Bim = 25.56, Bimax = 35.63, and the standard deviation σn-1 = The (tcalc) with or without the use of dT(1)/dt at x = 0 when using model (2) are similar or equal. Therefore, it is generally unnecessary to use dT(1)/dt at x = 0.
The (tcalc) from model (2) using dT(n)/dt at x = L is generally different from the value calculated without this derivative. Therefore, it is important to use dT(n)/dt at x = L
AN US
in model (2).
Generally, the oscillations are smaller with k((1/2)(Tn-1+Tn-2)) at x = L in model (2), than those with k(Tn-2).
The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
M
commercial calculation of (tcalc) agrees well with the experimental data (texper).
ED
4.2. Simulations using experimental data from Bazan and Mascheroni
The highlights of the prediction of (tcalc) using data from Bazan and Mascheroni (1984)
PT
by employing case (B) are as follows: Table 2 shows the results of the calculations with model (2) using k(Tn-2) at x = L.
CE
The Bif obtained were Bimin = 1.55, Bim = 1.96, Bimax = 2.35, and σn-1 = 0.33. The typical runtime is 3 s.
AC
The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or commercial calculation of (tcalc) agrees well with the experimental data (texper).
4.3. Simulations using experimental data from Cleland and Earle
The highlights of the prediction of (tcalc) using the data reported by Cleland and Earle (1977a; p. 1392) by employing case (B) are as follows:
ACCEPTED MANUSCRIPT The Bif obtained for tylose gel slabs were Bimin = 0.22, Bim = 3.56, Bimax = 20.78, and σn-1 = 5.69. The typical runtime is 4 s. The Bif obtained for mashed potato slabs were Bimin = 0.22, Bim = 1.81, Bimax = 4.21, and σn-1 = 1.59. The typical runtime is 3 s. The Bif obtained for minced lean beef slabs were Bimin = 0.27, Bim = 2.83, Bimax = 10.29, and σn-1 = 3.80. The typical runtime is 3 s. The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
CR IP T
commercial calculation of (tcalc) agrees well with the experimental data (texper).
4.4. Simulations using experimental data from Hense and Kieckbusch
Kieckbusch (1991) are the following:
AN US
The relevant aspects of the prediction of (tcalc) using the data from Hense and
The Bif obtained were Bimin = 4.05, Bim = 5.10, Bimax = 6.28, and σn-1 = 1.03. The typical runtime is 3 s.
The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
M
commercial calculation of (tcalc) agrees well with the experimental data (texper).
ED
4.5. Simulations using experimental data from Hung and Thompson
The relevant aspects of the prediction of (tcalc) using the data from Hung and Thompson
PT
(1983) are the following:
Only one experiment on tylose gel slabs, with Bif = 0.09 and a water content of
CE
0.751, required a large n and the runtime in case (B), using the data from Hung and Thompson (1983).
AC
The Bif obtained for tylose gel slabs were Bimin = 0.08, Bim = 1.15, Bimax = 3.37, and σn-1 = 0.99. The typical runtime is 5 s.
The Bif obtained for lean beef slabs were Bimin = 0.08, Bim = 1.28, Bimax = 3.30, and σn-1 = 1.09. The typical runtime is 3 s. The Bif obtained for ground beef slabs were Bimin = 0.06, Bim = 1.17, Bimax = 3.10, and σn-1 = 1.04. The typical runtime is 3 s. The Bif obtained for carp slabs were Bimin = 0.06, Bim = 1.16, Bimax = 3.15, and σn-1 = 1.08. The typical runtime is 3 s.
ACCEPTED MANUSCRIPT The Bif obtained for mashed potato slabs were Bimin = 0.06, Bim = 1.11, Bimax = 3.03, and σn-1 = 1.03. The typical runtime is 3 s. The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or commercial calculation of (tcalc) agrees well with the experimental data (texper).
4.6. Simulations using experimental data from Lescano
CR IP T
The relevant aspects of the prediction of (tcalc) using the data from Lescano (1973) are the following:
The Bif obtained were Bimin = 0.86, Bim = 3.06, Bimax = 6.18, and σn-1 = 1.59. The typical runtime is 3 s.
The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
AN US
commercial calculation of (tcalc) agrees well with the experimental data (texper).
4.7. Simulations using experimental data from Mascheroni and Calvelo
Calvelo (1982) are the following:
M
The relevant aspects of the prediction of (tcalc) using the data from Mascheroni and
ED
The Bif obtained for lean beef slabs were Bimin = 2.82, Bim = 5.62, Bimax = 8.76, and σn-1 = 1.65. The typical runtime is 3 s. The Bif obtained for lean beef slabs with heat transfer parallel to the fibers were Bimin
PT
= 1.39, Bim = 4.70, Bimax = 7.36, and σn-1 = 2.66. The typical runtime is 3 s. The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
CE
commercial calculation of (tcalc) agrees well with the experimental data (texper).
AC
4.8. Simulations using experimental data from Michelis and Calvelo
The relevant aspects of the prediction of (tcalc) using the data from Michelis and Calvelo (1983) are the following: The Bif obtained were Bimin = 2.81, Bim = 5.18, Bimax = 7.36, and σn-1 = 1.97. The typical runtime is 3 s. The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or commercial calculation of (tcalc) agrees well with the experimental data (texper).
ACCEPTED MANUSCRIPT
4.9. Simulations using experimental data from Pham and Willix
The highlights for the prediction of (tcalc) using the data from Pham and Willix (1990) are as follows: The Bif obtained were Bimin = 1.17, Bim = 7.57, Bimax = 31.93, and σn-1 = 6.59. The typical runtime is 4 s.
CR IP T
The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
commercial calculation of (tcalc) agrees well with the experimental data (texper).
5. Conclusions
AN US
A numerical model has been developed for predicting the freezing process of foods or analogues. The finite difference method (FDM) was employed to calculate the temperature distribution in a food T(x,t). The freezing times predicted (tcalc) with model (2) using the properties of tylose gel were compared with published experimental data
M
(texper).
In summary, the main conclusions of this study are as follows:
ED
It is important to use the derivative dT(n)/dt at x = L in model (2) in relation to the use of dT(1)/dt at x = 0 in the same model. The use of k(Tn-2) at x = L with model (2) and case (B) in the industrial or
PT
commercial calculation of (tcalc) is adequate. Simulations using k(Tn-2) at x = L with model (2) and case (B) were performed
CE
obtaining (tcalc) and compared with 227 experimental data points (texper), exhibiting good agreement. The parameters obtained for 227 experimental data points were Emin
AC
= -6.97 %, Emean = 0.06 %, Emax = 8.08 %, and n-1 = 2.19 %.
REFERENCES
Albasiny, E.L., 1956. The solution of non-linear heat-conduction problems on the Pilot Ace. Proceedings Inst. Electric Eng. 103, Part B, issue 1, 158-162. Bazan, H.C., Mascheroni, R.H., 1984. Heat transfer with simultaneous change of phase in freezing boned mutton. Lat. Am. J. Heat Mass Transfer. 8, 55-76.
ACCEPTED MANUSCRIPT Bonacina, C., Comini, G., 1973. On the solution of nonlinear heat conduction equations by numerical methods. Int. J. Heat Mass Transfer. 16, 581-589. Carslaw, H.S., Jaeger, J.C., 1959. Conduction of Heat in Solids, Clarendon Press. Cleland, A.C., 1977. Heat transfer during freezing of foods and prediction of freezing times. Ph. D. thesis, Massey University, New Zealand. Cleland, A.C., Earle, R.L., 1977a. A comparison of analytical and numerical methods of predicting the freezing times of foods. J. Food Sci. 42(5), 1390-1395.
CR IP T
Cleland, A.C., Earle, R.L., 1977b. The third kind of boundary condition in numerical freezing calculations. Int. J. Heat Mass Transfer. 20, 1029-1034.
Cleland, A.C., Earle, R.L., 1979. A comparison of methods for predicting the freezing times of cylindrical and spherical foodstuffs. J. Food Sci. 44, 958-963, 970.
Cleland, A.C., Earle, R.L., 1984. Assessment of freezing time prediction methods. J.
AN US
Food Sci. 49, 1034-1042.
Cleland, D.J., Cleland, A.C., Jones, R.S., 1994. Collection of accurate experimental data for testing the performance of simple method for food freezing time prediction. J. Food Proc. Eng. 17, 93-119.
M
Crank, J., Nicolson, P., 1947. A practical method for the numerical integration of solutions of partial differential equations of heat conduction type, Proceedings Cambridge Philosophical Society, 43, 50–67.
ED
Creed, P.G., James, S.J., 1983. The freezing times of liver in a vertical plate freezer. Progress in Refrigeration Science and Technology, Proceedings of the 16th
PT
International Congress of Refrigeration, Tome III, Paris, France, 4, 145-151. Creed, P.G., James, S.J., 1985. Heat transfer during the freezing of liver in a plate
CE
freezer. J. Food Sci. 50, 285-288, 294. Cutlip, M.B., Shacham, M., 1999. Problem Solving in Chemical Engineering with
AC
Numerical Methods, Prentice Hall. Dehghan, M., 2006. A computational study of the one-dimensional parabolic equation subject to nonclassical boundary specifications. Num. Meth. Partial Dif. Eq. 22, 220257.
Dima, J.B., Santos, M.V., Baron, P.J., Califano, A., Zaritzky, N.E., 2014. Experimental study and numerical modeling of the freezing process of marine products. Food Bioproducts Processing. 92, 54-66. Ferreira, S.R., Rojas, L.O.A., Souza, D.F.S., Oliveira, J.A., 2016. Freezing time of an infinite cylinder and sphere using the method of lines. Int. J. Refrigeration. 68, 37-49.
ACCEPTED MANUSCRIPT Ferreira, S.R., 2017. Freezing time of a slab using the method of lines. Int. J. Refrigeration. 75, 77-94. Hense, H., Kieckbusch, T.G., 1991. Freezing of shark: I Experimental results. In Proceedings of the IV Latin American Congress of Heat and Mass Transfer. 112 Hung, Y.C., Thompson, D.R., 1983. Freezing time prediction for slab shape foodstuff by an improved analytical method. J. Food Sci. 48, 555-560.
Corporate Headquarters, Houston, Texas.
CR IP T
Imsl, 2006. Fortran numerical library user‟s guide, version 6.0, Visual Numerics
Lees, M., 1966. A linear three level difference scheme for quasi-parabolic equations. Maths. Comput. 20, 516–522.
Lescano, C.E., 1973. Predicting freezing curves in codfish fillets using the ideal binary solution assumption. MS thesis, Michigan State University.
AN US
Luikov, A.V., 1968. Analytical Heat Diffusion Theory, Academic Press.
Mannapperuma, J.D., Singh, R.P., 1988. Prediction o freezing and thawing times of foods using a numerical method based on enthalpy formulation. J. Food Sci. 53, 626630.
M
Mannapperuma, J.D., Singh, R.P., 1989. A computer-aided method for prediction of properties and freezing/thawing times of foods. J. Food Eng. 9, 275-304. Mascheroni, R.H., Calvelo, A., 1982. A simplified model for freezing time calculations
ED
in foods. J. Food Sci. 47, 1201-1207. Michelis, A., Calvelo, A., 1983. Freezing time predictions for brick and cylindrical-
PT
shaped foods.J. Food Sci. 48, 909-913, 934. Northrop, P.W.C., Ramachandran, P.A., Schiesser, W.E., Subramanian, V.R., 2013. A
CE
robust false transient method of lines for elliptic partial differential equations. Chem. Eng. Sci. 90, 32-39.
AC
Pham, Q.T., 1985. A fast, unconditionally stable finite-difference scheme for heat conduction with phase change. Int. J. Heat Mass Transfer. 28(11), 2079-2084.
Pham, Q.T., 2008. Modelling of freezing processes, in: Evans, J.A. (Ed.), Frozen Food Science and Technology. Blackwell, Oxford, U.K., pp. 51-80.
Pham, Q.T., 2014. Food Freezing and Thawing Calculations, Springer, New York. Pham, Q.T., Willix, J., 1990. Effect of Biot number and freezing rate on accuracy of some food freezing time prediction methods. J. Food Sci. 55(5), 1429-1434. Sadiku, M.N.O., Obiozor, C.N., 2000. A simple introduction of the method of lines. Int. J. Electrical Eng. Educ. 37(3), 282-296.
ACCEPTED MANUSCRIPT Salvadori, V.O., 1994. Heat transfer during freezing, storage and thawing of foods. Ph. D. thesis, University National of La Plata, Argentina. (In Spanish.) Schiesser, W.E., 1991. The Numerical Method of Lines: Integration of Partial Differential Equations, Academic Press, San Diego. Schiesser, W.E., Silebi, C.A., 2009. Computational Transport Phenomena: Numerical Methods for the Solution of Transport Problems, Cambridge University Press. Schwartzberg, H.G., Singh, R.P., Sarkar, A., 2007. Freezing and thawing of foods –
CR IP T
computation methods and thermal properties correlation, in: Yanniotis, S., Sunden, B. (Eds.), Heat Transfer in Food Processing: Recent Developments and Applications, WIT Press, pp. 61-99.
Singh, R.P., Heldman, D.R., 2009. Introduction to Food Engineering, Elsevier.
Succar, J., Hayakawa, K., 1983. Empirical formulae for predicting thermal physical
AN US
properties of food at freezing or defrosting temperatures. Lebensm.-Wiss.u-Technol. 16(6), 326-331.
Wilson, H.A., Singh, R.P., 1987. Numerical simulation of individual quick freezing
AC
CE
PT
ED
M
spherical foods. Int. J. Refrigeration. 10, 149-155.
ACCEPTED MANUSCRIPT
i =1 T(1,t)
2 T(2,t)
3 T(3,t)
(i-1) T(i-1,t)
i T(i,t)
(i+1) T(i+1,t)
CR IP T
Fig. 1 – Grid of the half-thickness of a slab.
Table 1 – Computer program and subroutine for model (2) using tylose gel properties by a selected experiment. Program SlabDerivativeCaseB
M
PT
ED
Write(10,10) tend, (y(i), i = 1, n) t = tend tend = t + step End Do 10 Format (
AN US
Include 'link_fnl_shared.h' Use IVMRK_INT Implicit none Integer i, j, ido, n, np Parameter (n = 7, np = 5000) real*8 t, tend, Y(n), Yprime(n), step External Fcn ! Set initial values. t = 0.d0 Y(1:n) = 10.5 ! Defining final time and integration step. tend = 10000 step = (tend -t)/dble(np) tend = step ! Creating an output data file. Open (10, file = "Results.txt", status = "unknown") Write(10,10) t, (y(i), i = 1, n) Ido = 1 ! Normally, the initial call is made with IDO = 1. Do j = 1, np Call DIvmrk(ido, n, fcn, t, tend, y, yprime)
CE
!Ym = T(average) selected for calculation in x = L and i = n. !Ym = Y(n-1) Ym = Y(n-2) !Ym = (1/2)*(Y(n-2) + Y(n-1)) !Ym = (1/3)*(Y(n-2) + Y(n-1) + Y(n))
AC
If (Y(1).GT.-0.6) Then k = 0.49 Rho = 1006.0 Cp = 3688.0 Else If (Y(1).LE.-0.6) Then k = 1.573 - 0.00773*Y(1) + 0.653/Y(1) Rho = 937.908 - 0.154*Y(1) - 40.80/Y(1) Cp = 1418 + 122173/((-Y(1))**1.696) End If ! Evaluate derivative of temperature in x = 0 and i = 1: !Y(1) = Y(2) Yprime(1) = k*(Y(2)- 2*Y(2) + Y(3))/(Rho*Cp*(dx**2)) Do i = 2, n - 1 If (Y(i).GT.-0.6) Then k = 0.49 Rho = 1006.0 Cp = 3688.0 Else If (Y(i).LE.-0.6) Then k = 1.573 - 0.00773*Y(i) + 0.653/Y(i) Rho = 937.908 - 0.154*Y(i) - 40.80/Y(i) Cp = 1418 + 122173/((-Y(i))**1.696) End If
i=n T(n,t)
ACCEPTED MANUSCRIPT
! Evaluate derivatives of temperature dT/dt: Yprime(i) = k*(Y(i-1) - 2*Y(i) + Y(i+1))/(Rho*Cp*(dx**2)) End Do If (Ym.GT.-0.6) Then k = 0.49 Rho = 1006.0 Cp = 3688.0 Else If (Ym.LE.-0.6) Then k = 1.573 - 0.00773*Ym + 0.653/Ym Rho = 937.908 - 0.154*Ym - 40.80/Ym Cp = 1418 + 122173/((-Ym)**1.696) End If
CR IP T
! Evaluate derivative of temperature in x = L and i = n: !k*(Y(n)-Y(n-1)) = dx*h*(Ya-Y(n)) Yprime(n)= (k*(Y(n-2)-Y(n-1)) + dx*h*(Ya-Y(n)))/(Rho*Cp*(dx**2.0)) End subroutine
for slabs using case (B) and model (2).
AN US
Table 2 – Summary of comparisons between experimental and predicted freezing times
Material and model
Evaluation of (k) in x = L
Tc (°C)
Runs
Xmin (%)
Xmean (%)
Xmax (%)
n-1 (%)
Bazan and Mascheroni (1984)
Boned mutton
k(Tn-2)
-18
4
-1.25
-0.21
0.22
0.70
Cleland and Earle (1977a)
Tylose gel
k(Tn-2)
-10
43
-3.44
0.89
4.98
2.08
Cleland and Earle (1977a)
Mashed potato
k(Tn-2)
-10
6
-2.05
0.00
3.21
1.78
Cleland and Earle (1977a)
Minced lean beef
k(Tn-2)
-10
6
-1.98
0.06
2.31
1.68
Creed and James (1983; 1985)
Pig liver
k(Tn-2)
-7
44
-3.87
0.00
5.14
1.91
Pig liver
k((1/2)(Tn-2+Tn-1))
-7
44
-2.94
-0.20
1.69
0.76
Shark meat
k(Tn-2)
-10
6
-3.44
-0.83
1.33
1.71
ED
PT
CE
Hense and Kieckbusch (1991)
M
Source of data
Tylose gel
k(Tn-2)
-18
23
-6.97
-0.79
2.24
2.00
Hung and Thompson (1983)
Lean beef
k(Tn-2)
-18
9
-2.09
0.34
2.23
1.24
Hung and Thompson (1983)
Ground beef
k(Tn-2)
-18
9
-4.35
-0.89
5.24
2.70
Hung and Thompson (1983)
Carp
k(Tn-2)
-18
9
-5.37
-0.59
2.01
2.29
Hung and Thompson (1983)
Mashed potato
k(Tn-2)
-18
9
-4.20
-0.37
2.22
2.18
Lescano (1973)
Codfish fillet
k(Tn-2)
-18
8
-5.65
0.42
2.34
2.72
Mascheroni and Calvelo (1982)
Lean beef
k(Tn-2)
-18
10
-3.11
0.22
3.34
2.38
AC
Hung and Thompson (1983)
ACCEPTED MANUSCRIPT
Lean beef (heat transfer parallel to the fibers)
k(Tn-2)
-18
4
-3.82
0.50
8.08
5.48
Michelis and Calvelo (1983)
Lean beef
k(Tn-2)
-18
5
-3.10
-0.13
3.34
2.98
Pham and Willix (1990)
Tylose gel
k(Tn-2)
-18
32
-3.81
0.20
5.53
2.28
AC
CE
PT
ED
M
AN US
CR IP T
Mascheroni and Calvelo (1982)