Optimal dynamic heat generation profiles for simultaneous estimation of thermal food properties using a hotwire probe: Computation, implementation and validation

Journal of Food Engineering 84 (2008) 297–306 www.elsevier.com/locate/jfoodeng

Optimal dynamic heat generation profiles for simultaneous estimation of thermal food properties using a hotwire probe: Computation, implementation and validation Nico Scheerlinck a,*, Nahor H. Berhane b, Carmen G. Moles c, Julio R. Banga d, Bart M. Nicolaı¨ a a

Division MeBioS – Mechatronics, Biosystems and Sensors, Biosystems Department, Katholieke Universiteit Leuven, Willem de Croylaan 42, B-3001 Leuven, Belgium b Department of Agricultural and Biological Engineering, Purdue University, 225 South University Street, West Lafayette, IN 47907-2093, USA c ABB Process Automation Division, BU Oil and Gas Strategic Research, Bergerveien 12, 1375 Asker, Norway d Process Engineering Group, Instituto Investigaciones Marinas, Consejo Superior de Investigaciones Cientı´ficas, C/Eduardo Cabello 6, 36208 Vigo, Spain Received 15 September 2006; received in revised form 8 May 2007; accepted 12 May 2007 Available online 24 May 2007

Abstract The methodology of dynamic optimal experimental design was applied to estimate the thermal conductivity and the volumetric heat capacity of conduction heated foods simultaneously. For this purpose, a hotwire probe was used to generate heat in the thermal center of a food product. The temperature response was measured at a specified location inside the food. In order to minimize the effect of experimental input errors on the estimates, the information content of the measured temperature profiles was maximized. Optimal (dynamic) heat generation profiles as well as an optimal location to measure the temperature response in the food product were computed based on modern global optimization algorithms. These heating profiles were implemented and it was shown that the methodology guarantees an optimal regime for simultaneous parameter estimation. Pulse heating with two peaks – one at the beginning and one at the end of the experiment – was found to be an optimal heating profile. The optimal measurement position was located at some distance from, and close to, the heating probe. The mathematical and numerical outcomes of the applied methodology were interpreted in a physical context. Validation experiments with a food simulator (Tylose) were done, and coefficients of variation of 4% and 2% were obtained for the thermal conductivity and volumetric heat capacity, respectively. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Hotwire probe; Optimal experimental design; Global optimization; Heating strategy

1. Introduction As a result of the new European directives in food treatment processes, the organizations of the sector call for more fundamental research on the techniques to apply and, more particularly, for optimal treatment policies for their products. The general goal is obvious, a higher knowledge and control of the food products and processes, both regarding microbiological safety and final quality. *

Corresponding author. Tel.: +32 16 32 05 90; fax: +32 16 32 29 55. E-mail address: [email protected] (N. Scheerlinck).

0260-8774/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.jfoodeng.2007.05.019

Since the final quality and marketing of food products are closely related to their thermal treatment history, knowledge of the thermal product properties is crucial for the design, modelling, simulation and evaluation of different thermal food process scenarios and heating strategies. Besides the estimation of thermal food properties based on chemical composition of the food using empirical equations (Miles, van Beek, & Veerkamp, 1983) or structured models (Carson, Lovatt, Tanner, & Cleland, 2006; Wang, Carson, North, & Cleland, 2006), various measurement methods exist involving steady-state and transient-state

298

N. Scheerlinck et al. / Journal of Food Engineering 84 (2008) 297–306

heat transfer methods. Among others the transient hotwire probe method is frequently used to measure the thermophysical properties of conduction heated foods, and emerged to be the method of choice because of its accuracy (Hooper & Chang, 1952). In this method, a needle with a heating wire is inserted in the food and the heating is started. The thermal conductivity can then be calculated from the time–temperature data measured in the probe or at some position in the food. Recent interests are directed towards the simultaneous estimation of the thermal conductivity and heat capacity by means of single or dual heat probe methods (Bristow, Bilskie, Kluitenberg, & Horton, 1995; Bristow, Kluitenberg, & Horton, 1994; Muramatsu, Tagawa, Kasai, Sakai, & Fukushima, 1999; Nahor, Scheerlinck, Verniest, De Baerdemaeker, & Nicolaı¨, 2001; Nahor, Scheerlinck, Van Impe, & Nicolaı¨, 2003; Tagawa, Kitamura, Muramatsu, & Murata, 1996). Both the quality and the quantity of the experimental data are commonly recognized as important properties, which contribute to the parameter estimation accuracy. High quality data are mostly acknowledged as experimental observations featuring low measurement error. Moreover, the persistency of excitation of the system – through heat generation – contributes to an accurate estimation of the involved parameters. This means that if the dissipated power by the hotwire probe fails to excite the food system, the measured signal may not be sufficiently rich to estimate the parameters simultaneously. Previous research showed that by using specific time varying heat generation profiles more accurate estimates were obtained as compared to the traditional method of constant heating (Nahor et al., 2003). Likewise, the degree of excitation may not be equal at every location in the food, leading to the existence of an optimal measurement position associated to a given heat generation profile applied (Nahor et al., 2003). Although evidence was shown that optimal dynamic heating generation profiles and optimal sensor locations must exist, it is still unclear how to optimize the information content of an experiment with respect to these factors and how to compute an optimal dynamic heat generation profile. In this contribution modern global optimization algorithms are used for this purpose and the resulting ‘‘theoretical” dynamic heating profiles are implemented in an experimental hotwire probe setup and evaluated on Tylose, i.e., a food analogue. In the discussion section, the mathematical and numerical outcomes of the methodology are interpreted in a physical context, and factors affecting the implementability of the heating profiles are addressed. 2. Materials and methods 2.1. Hot wire probe setup In order to estimate the thermophysical properties of conduction heated food simultaneously, a dual probe hot wire setup was designed and constructed (Fig. 1).

Fig. 1. (a) Schematic representation of hot wire probe set up and (b) details of the heating probe (dimensions in mm).

The hot wire probe was made from a stainless steel hypodermic needle 21G  200 (Terumo, The Netherlands) which encases an insulated constantan (Cu55/Ni45) resistance wire (Good fellow Cambridge Ltd., UK) with a diameter of 0.07 mm. Constantan was selected for its distinct characteristics that the resistance is practically independent of temperature. The remaining space in the needle hole was filled with a heat tolerant resin (810 electrical resin, Schotchcast, Germany). The probe was inserted at the geometrical centre of a container (D = 35 mm and L = 90 mm) filled with a food simulator medium. A calibrated type-T TeflonÓ-insulated thermocouple – with a diameter of 0.2 mm (Omega Engineering, Stanford, USA) and connected to an HP 34970A data acquisition system (Agilent technologies, USA) – was used to measure the temperature in the medium at a known distance from the heating probe. Electrical current was supplied to the heating probe by a DC power source (HP E3632A, Agilent technologies, USA), which is equipped with a multi-meter to measure the voltage and current. For the acquisition of the time temperature data and the applied power a data logger (HP 34970A, Agilent technologies, USA) was used. In order to control the hardware for the dynamic heat generation profile and the synchronization of the data acquisition from the different measuring devices, a program was written using the LabView (National instruments Co., Austine Texas) data acquisition environment. Tylose – a methylcellulose gel – was used as food simulator medium for the fact that its thermophysical properties are similar to most common food substances. Tylose was prepared by mixing Tylose powder with water to give final moisture content of 77% wet basis. The thermophysical properties of Tylose food simulant are well known, i.e., a thermal conductivity of 0.56 W m1 °C1 and a volumetric heat capacity of 3.712  106 W m3 (Cleland, 1990; Incropera & De Witt, 1990; Reidel, 1960). These literature values were used in the design methodology.

N. Scheerlinck et al. / Journal of Food Engineering 84 (2008) 297–306

2.2. Heat transfer model

CHAMPACK routines (Scheerlinck, Verboven, Fikiin, De Baerdemaeker, & Nicolaı¨, 2001; Scheerlinck et al., 2004).

Heat conduction in the hotwire probe setup was modelled by means of the Fourier equation (Bird, Stewart, & Lightfoot, 2002) oT ¼ kr2 T þ QðtÞ ot T ðt ¼ t0 ; x; y; zÞ ¼ T 0 ðx; y; zÞ qcðT Þ

ð1Þ ð2Þ

where qc is the volumetric heat capacity [J m3 °C1], T the temperature [°C], t the time [s], k the thermal conductivity [W m1 °C1] and Q(t) the time-dependent volumetric heat generation [W m3]. The initial condition for the Fourier equation (Eq. (2)) is described as a spatial function at t = 0. At the boundary, C, of the food container convection heat transfer to the environment was considered (Eq. (3)) k

oT ¼ hðT  T 1 Þ on

299

ð3Þ

where n is the outward normal to the container surface, h is the surface heat transfer coefficient [W m2 °C1] and T1 the temperature of the environment, i.e., room temperature [°C]. A schematic representation of the geometric model for the hot wire probe set up is given in Fig. 2. The detailed components of the heating probe, namely the heating wire, the filler resin, and the needle wall, were incorporated in the spatial domain making up the heat transfer problem. A 1-D finite element model consisting of 81 nodes and 80 1D linear iso-parametric elements was built for the numerical solution of the heat transfer equation (Eqs. (1) and (2)). For numerical details on the finite element method the reader is referred to, e.g., Zienkiewicz and Taylor (1994). The finite element model and algorithms to solve the governing systems of differential equations were implemented in MATLAB (The MathWorks Inc., Natick, Massachusetts) using a multistep ODE solver (Shampine & Reichelt, 1997) and

2.3. Optimal experimental design (OED) The aim of OED is to find optimal experimental and operational settings to estimate the thermophysical properties from the time–temperature history measured at a specific location in the food substance, at some distance from the hotwire probe. The focus here is to achieve an informative experiment by using an optimal heat generation profile, Q(t), and an optimal measurement location, where optimal here means maximal information content. In practice, identifiability of thermophysical parameters requires heat generation input profiles ensuring that the effect of measurement uncertainties are minimized, resulting in robust parameter estimates. In previous research (Nahor et al., 2001, 2003) it was shown that in case of traditional heating strategies with constant heat generation profiles and noise corrupted temperature responses an estimation procedure might result in parameter estimates which deviate from the real physical parameter values to a large extent. Aiming at a simultaneous estimation of conduction heat transfer parameters we are left with the question whether the numerical values of the estimates correspond to physically meaningful heat transfer parameters or not. In what follows we will apply optimal experimental design techniques to estimate the thermophysical properties of food substances with the objective of defining the experimental procedure such that the sensitivity of input errors is minimized and the parameters are estimated with greatest statistical confidence. Hereto, the entire thermal system must be excited in such a way that the information content of the measured temperature response is maximized (optimal). Optimization of experimental procedures include initial conditions, experimental inputs, location of the hotwire probe and the like. To estimate the thermophysical parameters an iterative procedure is used in such a way that the final parameter estimates minimize some identification cost functional. The cost or identification functional, J, quantifies the deviation between the model predictions and the experimental data (Eq. (4)) Z tf T J¼ ðy  ym Þ Qðy  ym Þ dt ð4Þ t0

Fig. 2. 1D axi-symmetric geometric model (dimensions in mm). Constantan wire ( ), resin ( ), stainless steel ( ), medium ( ) and copper ( ).

with y the vector of predicted temperatures at n measurement positions, ym the vector of measured temperatures at n measurement positions, Q a user supplied square weighing matrix, t0 and tf the initial and final time, respectively. Considering a perfect model structure, noise free experimental data, and an appropriate heating strategy, the parameter estimation procedure yields the correct and/or assumed (nominal) parameters, i.e., the thermophysical properties of the food test substance Tylose. However, depending on the magnitude of the noise involved and other experimental errors, a change in one parameter can

300

N. Scheerlinck et al. / Journal of Food Engineering 84 (2008) 297–306

sometimes be compensated almost completely by a proportional shift in another. The effect of experimental errors on the estimation of the parameters was illustrated in previous research (Nahor et al., 2001, 2003). Furthermore, several combinations of numerical values for the thermal conductivity and the heat capacity might result in the same numerical value for the thermal diffusivity, which is the main unknown in the heat conduction model and is also the characteristic parameter of diffusion type models. So, in order to obtain experimental data that contain sufficiently rich information to estimate the thermophysical parameters, the input signals, i.e., heat generation profiles, need to be designed in such a way that the information content of the experiment is optimal. In our approach we use a conceptual measure from information theory, i.e., the Fisher information matrix, F, to define the information content of an experiment (Lutwak, Yang, & Zhang, 2005; Munack, 1989): Z t f  T   oy oy F¼ Q dt ð5Þ op op t0 The matrix oy/op contains the sensitivity functions (change of model predicted temperatures with respect to the thermophysical parameters) and quantifies the dependence of temperature predictions on the parameter values relative to the nominal parameter set p*. The matrix, Q, is the inverse of the measurement error covariance matrix and, thus, reduces the impact of large noise corrupted measurements on the information content. The parameter vector p has two elements: the thermal conductivity and the volumetric heat capacity, p = [k, qc]. For the quantification of the quality of the parameter estimates usually scalar quantities of the Fisher information matrix, F, are used to reflect specific optimal design criteria. Several optimal design criteria have been suggested in previous research (Mehra, 1974; Munack, 1989; Walter & Pronzato, 1990). In this research the modified E-criterion was used, which minimizes the ratio of the largest (kmax) to the smallest (kmin) eigenvalue of the Fisher information matrix, F. The advantage of the modified E-criterion is that the minimum value is one. For the case of two parameters, the contours of the identification functional, J, are circles. The computed optimal heat generation profile, Q(t), will be optimal in the sense that the information content, F, of the experiment is maximized, i.e., the E-criterion is minimized. The sensitivity functions, (Eq. (5)), were calculated using a numerical scheme as described by Nahor et al. (2001). Algorithms to compute the Fisher information matrix and the E-criterion – based on the heat transfer model, (Eqs. (1)–(3)) – were developed and implemented using MATLAB (The MathWorks Inc., Natick, Massachusetts). 2.4. OED via global optimization Mathematically, the OED problem we are dealing with can be formulated as a dynamic optimization (optimal control) problem. Through a model-based approach we try to

find a time varying (dynamic) heat generation profile, Q(t), in such a way that the information content, F, of an experiment is maximal. The objective is to find a set of timevarying input variables (controls) for the experiments optimizing the quality of the estimated parameters in a statistical sense. This optimal control problem is subject to a set of equality constraints (the Differential-Algebraic Equations (DAEs) describing the system dynamics) and a set of inequality path and/or point constraints on the state variables, i.e., the temperature at a specific location in the food substance. These inequality constraints can be associated with restrictions concerning issues like practical implementation, safety and/or model validity. Here, we have considered the modified E-criterion as the performance index to minimize, using the heating profile, Q(t), as the control variable. Numerical solutions to this dynamic optimization problem can be obtained using direct methods, which transform the original problem into a non-linear programming (NLP) problem via parameterizations. However, due to the frequent non-smoothness of the cost functions, the use of gradient-based methods to solve this NLP might lead to local solutions (Banga, Versyck, & Van Impe, 2002). These authors have also shown, for the case of systems described by ordinary differential equations, that stochastic methods of global optimization can be used as robust alternatives. In this work, a general approach, based on the framework of control vector parameterization (Banga, BalsaCanto, Moles, & Alonso, 2003), was pursued by defining a general block-type heating profile through characteristic parameters such as the amplitude, duration and frequency of the peaks. This was accomplished by piecewise discretization of the input profile. Hereto, simple basis functions such as steps were introduced. This leads to a substantial number of degrees of freedom transforming the problem into a multi-dimensional non-linear optimization problem. For a given type of heat generation profile, Q(t), factors that influence the information content of an experiment are the duration and magnitude of the heat generation power. Moreover, as the degree of excitation might not be the same at all the points in the sample, the location at which the time–temperature measurement should be taken was also considered as an additional degree of freedom. This requires computationally more powerful algorithms to solve the parametric optimization problem under several imposed constraints, which may ensure global optimal solution. In the present work, two stochastic strategies for global optimization were employed, namely, the Integrated Controlled Random search method (Banga & Casares, 1987) and the Stochastic Ranking Evolution Strategy (Runarsson & Yao, 2000). Parametrization of the input heating profile was realized by discretizing the input heat generation profile into 10 steps. Only steps were considered because previous research revealed that general block-type profiles give more informative experiments (Nahor et al., 2003).

N. Scheerlinck et al. / Journal of Food Engineering 84 (2008) 297–306

Two cases were considered with respect to the measurement position, in the first case the measurement position was part of the parameter set to be optimized. Alternatively, the position was fixed to 4 mm in the second case. In the computations a uniform initial temperature for Tylose, T0, of 20 °C was considered. The ambient (room) temperature, T1, was set to 25 °C. The surface heat transfer coefficient, h, was set to 6.5 W m2 °C1; this value was estimated based on temperature measurements in the container filled with pure water at 20 °C. However, in hotwire probe experiments the surface heat transfer at room temperature does not significantly affect the temperature response in the sensor location. For the computation of the optimal heat generation profiles, the total duration of a heating experiment, tf, was fixed to 300 s. Additional constraints were applied on the magnitude of the heat generation in such a way that the temperature increase was restricted to some extend to avoid excessive overheating of the product. The optimal heat generation profiles, obtained from the solution of the global optimization problem, were implemented in the software, controlling the heat generation in the hotwire probe. Three – independent – experiments were done to estimate the thermophysical properties of Tylose samples. For each experiment, new Tylose samples were prepared. The hotwire probe and the thermocouple were positioned prior to the hardening of the sample. The hotwire probe was positioned in the center of the sample (Fig. 1). The thermocouple, used for the measurement of the temperature response, was placed in its optimal position, which was obtained from the solution of the global optimization problem. Prior to experimentation, the Tylose samples were conditioned at constant room temperature (T0 = T1 = 25 °C). When the heating experiment was finished, the exact position of the thermocouple was determined with a digital calliper (Mitutoyo Ltd., UK, ±0.01 mm). The estimated thermophysical parameters were compared with literature values (Cleland, 1990; Incropera & De Witt, 1990; Reidel, 1960). An error analysis, based on simulation and experimental results, was performed to check the accuracy and robustness of the developed methodology. 3. Results The optimal heat generation profiles obtained for the two cases are presented in Fig. 3. In the first case (a), where the optimization procedure was left to determine the optimal heat generation profile, Q(t), as well as the optimal measurement position, the a priori known global minimum has been attained, i.e., the modified E-criterion is equal to one. Since the modified E-criterion represents the ratio of the largest to the smallest eigenvalues of the Fisher Information matrix, this ratio can not have a value smaller than one. And, therefore, the information content of an experiment for which the modi-

301

Fig. 3. Optimal heating profile from parameterized optimization: (a) DOF = 11 – measurement location = 2.25 mm (free) – modified Ecost = 1.0; (b) DOF = 10 – measurement location = 4.0 mm (fixed) – modified E-cost = 1.32. Pulse heating profiles (left scale) and temperature response (right scale).

fied E-criterion is equal to one, is not only optimal, but maximal in the full sense. When this global minimum is reached, the associated heat generation profile and measurement position are called globally optimal or optimal in the full sense. Pulse heating with two peaks, one at the beginning and one at the end of the experiment, was identified as a globally optimal profile with the optimal sensor position located at 2.25 mm from the centre of the probe. In the second case (b), where the measurement position is imposed (fixed) at 4.0 mm from the center of the hot wire probe, the optimal solution was found at a slightly higher value than the global minimum. A single heat pulse at the beginning followed by a ‘‘gradual” decrease of heat generated was found as an optimal profile. When the measurement position is (virtually) altered, a different heating profile is obtained. In Table 1 the results of an error analysis are shown. Three types of errors were considered acting on the model inputs:

302

N. Scheerlinck et al. / Journal of Food Engineering 84 (2008) 297–306

Table 1 Error analysis, based on simulations with an optimal heating profile and optimal measurement position (Fig. 3a) Characterization of errors

Effect on estimates of thermal conductivity, k [W/m/°C]

Effect on estimates of volumetric heat capacity, qc [MJ/m3/°C]

Systematic error on the temperature response T(t) ± 1 °C 0.8% (relative error)

0.2% (relative error)

Noise on temperature response (100 repetitions) T(t) ± N(0, 0.03)a k ± 0.03 (95% confidence interval)

qc ± 0.19 (95% confidence interval)

Error on sensor position readingsb Exact 0.7% position  0.01 mm Exact 0.7% (relative error) position + 0.01 mm

1.6% 0.6% (relative error)

a

Noise generated from a gaussian distribution N(l, r), where r was determined from temperature measurements in temperature conditioned cool rooms. b Error bound (accuracy) of calliper = 0.01 mm.

(i) Systematic errors on the temperature response, (ii) Noise on the temperature response, and (iii) Errors related to the measurement of the sensor position, in which the temperature response is recorded. The effect of these errors on the parameter estimates was assessed through virtual heating experiments, simulated on a computer, and using the optimal input settings as specified by the outcome of the global optimization algorithm (Fig. 3a). Errors related to the heat transfer model, and its geometric representation, are not considered here, for it was shown in previous work that the model predicts well the actual temperature response, within the error limits of temperature measurements (Nahor et al., 2003). The lower and upper error limits of the parameter estimates are as follows: (i) A systematic error of 1 °C on the temperature response, i.e., T(t) ± 1 °C, causes relative errors of 0.8% and 0.2% on the estimated values for the thermal conductivity and volumetric heat capacity, respectively. Since the error on the temperature measurements is ±0.3 °C (this value was obtained from the calibration of the thermocouple), the relative errors are overestimated, and represent a worse case scenario. (ii) The effect of noise on the temperature response results in error bounds (at a 95% confidence level) of k ± 0.03 (W m1 °C1) and qc ± 0.19 3 1 (MJ m °C ). The 95% confidence intervals on the estimates were computed based on 100 simulated heating experiments, in which noise was added to the temperature response. The noise was generated from a gaussian distribution with zero mean and a standard deviation of 0.03 °C, which was identified from measurements in temperature conditioned cool

rooms. However, given the perfect fitting between sensor and material, noise on the temperature response measured in the samples is significantly filtered out due to the diffusive nature of conduction heat transfer. (iii) An error of 0.01 mm on the reading of the sensor position causes a 0.7% relative error on the estimated value for the thermal conductivity, and lower and upper bounds of 1.6% and 0.6% (relative errors) for the volumetric heat capacity. In order to validate the optimal experimental design and global optimization methodologies, the heating profile with global minimum of unity (modified E-criterion = 1) was implemented. An exact implementation of the ‘‘computed” optimal heating profile turned out to be difficult because of practical limitations related to the control of the hardware involved and the acquisition frequency in the experimental set up. For this reason, the theoretically obtained time intervals were adjusted in such a way that the generated heating profile was close to the optimal one (Fig. 4). A small deviation between the actual and optimal sensor position has no influence on the estimation procedure of the parameters, since the exact sensor location was measured when heating experiments were finished, with an accuracy of 0.01 mm. A practical implementation of the theoretically optimal heat generation profile results in a different value for the modified E-criteria (2.71 instead of unity). But, the methodology that we used here ensures that the effect of experimental errors on the parameter estimates is minimal (Fig. 5). This is reflected in the cone-like structure of the functional, J, which predicts only one combination of parameters, irrespective of the initial guess in the parameter estimation procedure.

Fig. 4. Implemented optimal heating profile. Pulse heating profile (left scale) and temperature response (right scale). Optimal measurement location: 2.25 mm; actual measurement location: 1.55 mm.

N. Scheerlinck et al. / Journal of Food Engineering 84 (2008) 297–306

303

The parameter estimates and their statistics – obtained from three replicate experiments, in which the heat diffusion was controlled by the implementation of an optimal dynamic heating profile (Fig. 4) – are summarized in Table 2. The mean values of the parameter estimates are 0.556 (W m1 °C1), for the thermal conductivity, and 4.033 (MJ m3 °C1), for the volumetric heat capacity, with coefficients of variation of 4.12% and 2.30%, respectively. The parameter estimates (Table 2) are within the error bounds, as predicted by means of a numerical error analysis (Table 1). The relative deviation of the estimated parameters with respect to the literature values is 0.8% for the thermal conductivity, and 8.7% for the volumetric heat capacity (Table 2). These relative errors are similar to the order of magnitude of estimates obtained in previous research (Nahor et al., 2003). 4. Discussion

Fig. 5. Contour plots of the functional [log (J)] for (a) the theoretical and (b) the implemented optimal heat generation profile.

Table 2 Parameter estimates for Tylose – and their statistics – using an optimal heating profile and optimal sensor location (Fig. 4)

Estimates Experiment 1 Experiment 2 Experiment 3 Mean value % difference from literature valuesa Standard deviation Standard error Coefficient of variation a

Thermal conductivity, k [W/m/°C]

Volumetric heat capacity, qc [MJ/m3/°C]

0.530 0.563 0.574

3.930 4.060 4.110

0.556 0.8%

4.033 8.7%

0.023 0.013 4.12%

0.093 0.054 2.30%

Cleland (1990), Incropera and De Witt (1990), Reidel (1960).

The thermal conductivity of a material is a theoretical concept, and is used to describe how fast thermal energy flows through the material. Basically, the thermal conductivity describes a steady-state situation, and quantifies the heat flux (W m2) through the material when a constant temperature gradient (°C m1) exists over the material. The measurement of this material property is difficult because it usually takes a lot of time to reach a steady-state situation and, moreover, the temperature gradient over the material should be kept constant. When dynamic (nonsteady-state) heating processes are considered, another thermal property starts playing a role, i.e., the volumetric heat capacity. This property quantifies how much energy (J) is required for changing the temperature of a certain volume (m3) by 1°C. Experiments using a hotwire probe, in which dynamic heating strategies are involved, perform measurements relatively fast, and much faster as compared to steady-state measurements. Given the initial and geometrical conditions of the material, and the specifications of the dynamic heating profile, such as power versus time, the change of temperature with time at a specific location in the material is determined by both the thermal conductivity and the volumetric heat capacity. In order to obtain good and reliable estimates we are left with the following questions: ‘‘What is an optimal dynamic heating strategy?”, ‘‘What makes it optimal?”, ‘‘How does it look like?”, and ‘‘Is there an optimal measurement location?”. We answered these questions using the methodology of optimal experimental design in combination with global optimization techniques. This methodology enables a theoretical evaluation of different experimental protocols, using a hotwire probe set up, to estimate the thermophysical properties of food substances simultaneously. According to this methodology, a heating strategy is more appropriate when the effect of experimental input errors on the parameter estimates is minimized. The heating strategy is optimal when the information content of

304

N. Scheerlinck et al. / Journal of Food Engineering 84 (2008) 297–306

an experiment is maximal. In order to increase the information content, the system needs dynamic inputs to excite the full dynamics of the system. In the problem at hand, the thermal system must be excited in such a way that the temperature response contains sufficient information in order to extract the real quantitative values of both thermal conductivity and the volumetric heat capacity. When the sensor, for measuring the temperature response, was placed at a fixed and relatively large distance from the hotwire probe (4 mm, see Fig. 3b), the optimal heating profile consists of three subsequent heat pulses with decreasing power. The heat pulses are situated in the beginning of the experiment, for it takes some time to generate a temperature response at the measurement location, from which the information content is sufficiently rich. In this case, the thermal system is optimally excited by a decreasing level of heat supply, addressing the sensitivity of the system with respect to its thermal properties in an optimal way, i.e., no better excitation of the system exists. Although the heating profile is optimal in the sense that it maximizes the information content of the experiment, the absolute global minimum of unity for the modified E-criterion is not reached (E-criterion = 1.32). Since the final time of the experiment was constrained (tf = 5 min) and the sensor location was fixed at a relatively large distance from the probe, no extra information could be gained from introducing an additional heat pulse in the middle or at the end of the experiment. This result suggests that a randomly chosen and fixed sensor location does not guarantee optimal heat generation profiles in the full sense. When the temperature response is measured in an optimal sensor location, and given an experimental time constraint of five minutes, the optimal heating profile consists of two non-subsequent heat pulses: one pulse at the beginning of the experiment, and one pulse at the end of the experiment. Moreover, the optimal location for measuring the temperature response is close to the hotwire probe (2.25 mm, see Fig. 3a). This explains the reason of having an additional heat pulse at the end of the experiment, instead of introducing all heat pulses at the beginning, for it takes now less time for the heat to flow towards the sensor location and to generate a temperature response. First, the thermal system is excited by one pulse of heat supply. Then, a time period follows in which the system develops a heating and cooling response, characterizing the system dynamics, and addressing the sensitivity of the system with respect to its thermal properties. Finally, an additional heat pulse is used to increase and fully maximize the information content of the experiment, especially with respect to the speed of heat flow, which is required to introduce a temperature response in the (optimal) sensor location. This result shows that there exist an optimal heating profile which generates a temperature response with maximal information content (modified E-criterion = 1), provided that the temperature

response is measured in an optimal sensor position. Moreover, the outcome of the optimization problem shows that this optimal sensor position exists. So, the optimal heating strategy is to apply a short pulse in the beginning of an experiment and to extract as much information as possible from the pulse response. In general, a short pulse response contains more information as a step response, i.e., a constant heating regime (Nahor et al., 2001, 2003). Towards the end of the experiment a second pulse is introduced to extract more information from the rate of consequent temperature rise. Since the second pulse is right at the end it has the same effect as a very short step, and little information is lost from the pulse response. Using a purely experimental approach, Bristow et al. (1994) have reported that for the short duration pulse heating, no statistically significant difference was observed in the estimates by changing the magnitude of the dissipated heat. At the contrary, they found out that by varying the pulse width a significant difference in estimates were obtained especially in the thermal conductivity. They reported that the reasons for this are unclear. In light of the results we obtained in this work and also from the results of our previous works (Nahor et al., 2001, 2003), this can be explained as follows. An increasing pulse width changes the heating profile from pulse-type to block-type heating. Further, if the pulse width is increased to the final time of the experiment, it results in a constant heating profile. For a given measurement position the information content of the heat generation profile decreases as we go from a short duration heat pulse to a constant heating profile. Moreover, the optimal measurement position is located at a distance from the centre of the heat source, which is smaller for a short duration pulse-type heating profile and larger for a constant heating profile. Based on these insights we can attribute the change in estimates, as observed by Bristow et al. (1994) for pulse-type heat generation profiles with increasing pulse width, to a decrease of the information content of the experiments using these pulse-type heat generation profiles. Bilskie, Horton, and Bristow (1998) report a maximum deviation of 5.1 and 5.9% for Al2O3 and glycerol, respectively using their dual probe method. These values are comparable to the accuracy that we obtained (Table 2). It should be noted that the assumed values for the thermal material properties, as used in our design methodology, are values that we obtained from the literature. These values were estimated using other methods than we use or they were calculated based on their chemical composition. So, the presupposition that these parameters are exact, and thus the only correct physical values, can be questioned. Nevertheless, the methodology, as outlined in this work, ensures that the effect of experimental input errors on the parameter estimates are minimal (Fig. 5), irrespective of the nature of the parameter estimation procedure or initial guess. A deviation from the global minimum of unity observed on the E-criterion of the implemented dynamic heating

N. Scheerlinck et al. / Journal of Food Engineering 84 (2008) 297–306

profile is attributed to the ‘implementability’ of the ‘‘computed” optimal heating profiles. This means that practical limitations of experimental set ups, with respect to the frequency of measurement and control of hardware, might lead to substantial loss of information due to the inability to capture the full dynamics of ‘‘theoretical” heat generation profiles, especially when they are changing fast in small time intervals. This leads to loss of accuracy in the information content of the experiment in spite of the fact that the given profile was proven to be optimal having – theoretically – maximum information content. In the application of the methodology to estimate the parameters for Tylose, there was no major problem encountered in the implementation of the dynamic heating profile, but, at the contrary, it was found very difficult to position the thermocouple at the exact optimal measurement location, which is the major factor contributing to the deviation of the actual global minimum from the design global minimum of unity. This deviation caused no errors in the parameter estimation procedure to validate the methodology, because the exact sensor position was measured a posteriori. In practical or out-of-laboratory situations – e.g., online measurement and estimation of thermal food properties – the use of special equipment is recommended, or should be developed, to position the sensors at the optimal measurement locations, and within the error bounds as shown in the error analysis (Table 1). As a proof of principle, in this work, application of optimal experimental design methodologies was demonstrated using a food substance simulator, i.e., Tylose. When dealing with real food substances non-homogeneity and anisotropy come into play. This might require the reconsideration of the assumptions used in building the model, which were based on the homogeneity of the food substance. To the authors knowledge, besides the use of suspensions of expanded polystyrene beads in a gel (Carson, Lovatt, Tanner, & Cleland, 2004), there is no standardized heterogeneous food simulant available. Moreover, another aspect that needs attention in the process of optimal experimental design when dealing with real food substances is the design robustness. Although a deviation from literature values was noticed, the methodology was shown to be robust for the estimation of the thermophysical parameters of Tylose. In contrast, Versyck (2000) and Versyck and Van Impe (2000) have shown that optimal experimental design based on incorrect assumed parameter values can lead to loss of accuracy in the estimates. This can be surmounted through an iterative design cycle where the estimates from the first iteration cycle are considered as the correct (or nominal parameters) in the proceeding iteration cycle. The progress of the iteration can be stopped when a user specified criterion is met, for instance, when the estimated parameter sets in two subsequent iterations are sufficiently close. Finally, in the implementation of theoretically obtained optimal modus operandi, a minimal effect of experimental input errors on the accuracy of parameter estimates is guar-

305

anteed, but the statistical confidence level of the estimated parameters should be established through repetitions. 5. Conclusions The application of optimal experimental design, novel as it is in thermal food processing, provides a basis to evaluate theoretically the information content of a given experimental protocol and facilitates the design of informative experiments. Global optimization turned out to be a powerful methodology to derive optimal heat generation profiles for a hot wire probe setup to estimate thermal food properties simultaneously. The Fisher information matrix was applied as a measure to evaluate the information content of a thermal experiment. Pulse heating with two peaks – one at the beginning and one at the end of the experiment – was found to be optimal in the full sense. The optimal measurement position was located at a distance from, and close to, the heating probe. The coefficients of variation, obtained from the implementation of an optimal experiment protocol, were 4.2% and 2.3%, for the thermal conductivity and the volumetric heat capacity, respectively. A proof of principle was shown: given an experimental hotwire probe set up and a fixed final experimentation time, there exist an optimal dynamic heating strategy and an optimal sensor location. Here optimal means that the temperature response contains maximum information to estimate the thermal properties of the samples. However, in practical situations substantial loss of information might occur due to the inability of the hardware to capture the full dynamics of theoretical heating profiles, especially for heat generation profiles with steep changes. Acknowledgements Author Nico Scheerlinck is Postdoctoral Fellow with the Flemish Fund for Scientific Research (F.W.O. – Vlaanderen). The authors gratefully acknowledge the F.W.O. for funding this research. Authors Banga and Moles acknowledge financial support from the Spanish Government (MEC AGL2004-05206-C02-01/ALI) and Xunta de Galicia (PGIDIT05PXIC40201PN). References Banga, J.R., Casares, J.J. (1987). ICRS: application to a wastewater treatment plant model. In ‘‘Process Optimization’’ (IChemE Symposium Series No. 100), pp. 183–192. Banga, J. R., Balsa-Canto, E., Moles, C. G., & Alonso, A. A. (2003). Improving food processing using modern optimization methods. Trends in Food Science and Technology, 14(4), 131–144. Banga, J. R., Versyck, K. J., & Van Impe, J. F. (2002). Computation of optimal identification experiments for nonlinear dynamic process models: An stochastic global optimization approach. Industrial and Engineering Chemistry Research, 41, 2425–2430.

306

N. Scheerlinck et al. / Journal of Food Engineering 84 (2008) 297–306

Bilskie, J. R., Horton, R., & Bristow, K. L. (1998). Test of dual-probe heat-pulse method for determining thermal properties of porous materials. Soil Science, 163(5), 346–355. Bird, R. B., Stewart, W. E., & Lightfoot, E. N. (2002). Transport Phenomena (second ed.). New York: Wiley. Bristow, K. L., Bilskie, J. R., Kluitenberg, G. J., & Horton, R. (1995). Comparison of techniques for extracting soil thermal properties from dual-probe heat pulse data. Soil Science, 160(1), 1–7. Bristow, K. L., Kluitenberg, G. J., & Horton, R. (1994). Measurement of soil thermal properties with a dual-probe heat-pulse technique. Soil Science Society of America Journal, 58, 1288–1294. Carson, J. K., Lovatt, S. J., Tanner, D. J., & Cleland, A. C. (2004). Experimental measurements of the effective thermal conductivity of a pseudo-porous food analogue over a range of porosities and mean pore sizes. Journal of Food Engineering, 63(1), 87–95. Carson, J. K., Lovatt, S. J., Tanner, D. J., & Cleland, A. C. (2006). Predicting the effective thermal conductivity of unfrozen, porous foods. Journal of Food Engineering, 75(3), 297–307. Cleland, A. C. (1990). Food refrigeration processes. Analysis, design and simulation. England: Elsevier Science publishers Ltd.. Hooper, F. C., & Chang, S. C. (1952). Development of the thermal conductivity probe. Heating, Piping, and Air Conditioning, 24(10), 125–129. Incropera, F. P., & De Witt, D. P. (1990). Fundamentals of Heat and Mass Transfer (3rd ed.). New York, Chichester, Brisbane, Toronto, Singapore: John Willy and Sons, Inc. Lutwak, E., Yang, D., & Zhang, G. (2005). Cramer-Rao and momententropy inequalities for Renyi entropy and generalized fisher information. IEEE Transactions on Information Theory, 51, 473–478. Mehra, R. (1974). Optimal input signal for parameter estimation in dynamic systems: Survey and New Results. IEEE Transactions on Automatic Control, 19(6), 753–768. Miles, C. A., van Beek, G., & Veerkamp, C. H. (1983). Calculation of thermophysical properties. In R. Jowitt, F. Escher, B. Hallstro¨m, H. F. Meffert, W. E. L. Spiess, & G. Vos (Eds.), Physical properties of foods. UK: Applied Science Publishers. Munack, A. (1989). Optimal feeding strategy for identification of monodtype models by fed-batch experiments. In N. M. Fish, R. I. Fox, & N. F. Thomhill (Eds.), Computer applications in fermentation technology, modelling and control of biotechnological process (pp. 195–204). Amsterdam: Elsevier. Muramatsu, Y., Tagawa, A., Kasai, T., Sakai, H., & Fukushima, M. (1999). Measurements of thermophysical properties of wheat flour by transient heat flow method using a probe. Journal of the Japanese Society for Food Science and Technology, 46(11), 719–724.

Nahor, H. B., Scheerlinck, N., Van Impe, J., & Nicolaı¨, B. M. (2003). Optimization of the temperature sensor position in a hot wire probe set up for estimation of the thermal properties of foods using optimal experimental design. Journal of Food Engineering, 57(1), 103–110. Nahor, H. B., Scheerlinck, N., Verniest, R., De Baerdemaeker, J., & Nicolaı¨, B. M. (2001). Optimal design of dynamic experiments for the parameter estimation of conduction heated foods. Journal of Food Engineering, 48, 109–119. Reidel, L. (1960). Eine prufsubstanz fu¨r gefrierversuche. Kaltetechnik, 12, 222–226. Runarsson, T. P., & Yao, X. (2000). Stochastic ranking for constrained evolutionary optimization. IEEE Transactions on Evolutionary Computation, 4, 284–294. Scheerlinck, N., Marquenie, D., Jancsok, P. T., Verboven, P., Moles, C. G., Banga, J. R., et al. (2004). A model-based approach to develop periodic thermal treatments for surface decontamination of strawberries. Postharvest Biology and Technology, 34(1), 39–52. Scheerlinck, N., Verboven, P., Fikiin, K. A., De Baerdemaeker, J., & Nicolaı¨, B. M. (2001). Finite element computation of unsteady phase change heat transfer during freezing or thawing of foods by using a combined enthalpy and Kirchhoff transform method. Transactions of the American Society for Agricultural Engineers, 44(2), 429–438. Shampine, L. F., & Reichelt, M. W. (1997). The MATLAB ODE Suite. SIAM Journal on Scientific Computing, 18, 1–22. Tagawa, A., Kitamura, Y., Muramatsu, Y., Fusakaz A., Murata, S. (1996). Thermophysical properties of Adzuki beans by transient heat flow method using probe. ASEA Annual International Meeting. Paper # 966003. Versyck, K.J. (2000). Dynamic input design for optimal estimation of kinetic parameters in bioprocess models. Doctoral dissertation, No. 436. Faculty of Bioscience Engineering, Katholieke Universiteit Leuven, Belgium. Versyck, K. J., & Van Impe, J. F. (2000). Optimal design of a closed loop controller for estimation of parameter couples of microbial growth kinetics. Chemical Engineering Communications, 180, 39–59. Walter, E., & Pronzato, L. (1990). Qualitative and quantitative experiment design for phenomenological methods: A survey. Automatica, 26(2), 195–213. Wang, J. F., Carson, J. K., North, M. F., & Cleland, D. J. (2006). A new approach to modelling the effective thermal conductivity of heterogeneous materials. International Journal of Heat and Mass Transfer, 49(17–18), 3075–3083. Zienkiewicz, O. C., & Taylor, R. L. (1994). The finite element method: Basic formulation and linear problems. Berkshire, UK: McGraw-Hill.