Simulation of a commercial indirect UHT heat treatment

Koen Grijspeerdt, Leen Mortier, Jan De Block & Roland Van Renterghem

Centre of Agricultural Research
Department of Animal Product Quality
Brusselsesteenweg 370
B-9090 Melle
E-mail: K.Grijspeerdt@clo.fgov.be

Introduction

A model was implemented for a commercial UHT treatment system. The fact that there are 2 nested heat regeneration cycles in the system leads to excessive complication. Therefore, some simplifications were introduced. Several reactions are included in the model. The time-temperature profile is adequately reproduced by the model. Fouling was characterized by the pressure drop change during a process batch, sterilization was described by the inactivation of Bacillus stearothermophilus spores. Simulations clearly pointed to the bottlenecks with respect to fouling. However, the influence on sterilization is only minor. The influence of the mass flow on the formation of lactulose and hydroxymethylfurfural (HMF) on the other hand is very pronounced.
Some optimization strategies were tested out. It seems possible to lower the amount of fouling by adapting the UHT-conditions: a lower temperature at a longer residence time still guarantees sufficient bacterial inactivation and a lower fouling. Prolonging the residence time in the first intermediate holder can further reduce the fouling. Taking lactulose and HMF-formation into account, this latter strategy is preferred.

The heat treatment system

A schematic lay-out of the heat treatment system is shown in Figure 1.

The system has a fairly complicated lay-out. It is designed to obtain maximum heat recuperation, incorporating 2 nested heat regeneration cycles. This is in principle a very advanced design, but there are some problems with the heat distribution among the system. The heat link between the preheating and primary cooling module is not able to cool down the latter sufficiently. This then leads to a too severe heating in the first preheating (heat exchanger 1). The overall result is that the temperature in the first stage of the heat treatment is higher than the designed target values and that fouling can occur in this stage. The reason for the insufficient cooling capacity of the inner heat regeneration loop could be due to radiation effects. Radiation phenomena are not included in the mathematical model, because they require detailed descriptions of the geometry, which are not available. Besides, it would be too far off the primary goals of the modelling exercise in the first place.

Another peculiarity about the system is the use of 7-pipes heat exchangers. The product to be heated is sent through 7 small pipes with relatively small diameter and the heating medium flows through the external coil. This leads to higher heat exchange coefficients for the same heat exchanger length, but at the cost of a higher sensitivity towards fouling due to the relatively small diameters. Referring to Figure 1, heat exchanger 2, 3, 4 and Cooling 1 use  7 pipe systems. Holding tubes are 1-pipe systems.

The model

The heat treatment system is modelled as train of individual systems which can be heat exchangers, homogenization steps or holding tubes. In each of these systems chemical and microbiological reactions can occur.

Chemical reactions

Fouling

The model used is essentially the one presented by de Jong (1996). The chemical reactions leading to fouling are assumed to follow a two-stage pathway. βlg (denoted as in the model) first loses its ternary structure and becomes unfolded (). The unfolded βlg then is converted to aggregated βlg (), considered to play no significant role in the fouling process. De Jong (1996) determined this aggregation step to be a second order reaction. U can also react with milk constituents () to form aggregated milk constituents () which then can be adsorbed to the heat exchanger wall causing fouling (). The reaction scheme is shown in equation 1.

The rates of disappearance and formation are given by the corresponding reaction rate equations (equation 2).

The temperature dependencies of the reaction rates are described by an Arrhenius type law. The Arrhenius plot of the aggregation of βlg shows two distinct regions where the pre-exponential factor and the activation energies are different so the aggregation reaction has to be split over 2 temperature intervals. The fouling process itself () was experimentally verified to be a reaction-rate controlled process, not limited by mass transfer (de Jong, 1996). The experimental data were best described by a reaction of order 1.2 (equation 3). The values of the Arrhenius parameters were determined by de Jong (1996) and are summarized in Table 1.

It has to be stressed that the model of de Jong only focusses on fouling caused by βlg. There is another type of fouling,which becomes increasingly important at higher temperatures,namely fouling caused by precipitation of minerals. The nature of mineral fouling is quite different form protein fouling, and the interaction between the 2 is not well known. There are not many modelling efforts yet in this area, mainly because the underlying phenomena are very complex. The temperatures reached in UHT-treatment, as the one in this report, will give rise to a certain degree of mineral fouling. It should therefore be kept in mind that the predicted fouling will be an underestimation.

Other components

Some other chemical reactions were added, although these were not taken into consideration for optimization purposes. The kinetic expressions were derived at the Laboratory of Food Technology of the KU Leuven or derived from industrial simulation data. Some may be not very relevant to the case studied here, but were included in the model for the sake of completeness.

Alkalic phosphatase

Lactoperoxidase

Lactulose (kinetic expressions derived from industrial data)

Hydroxymethylfurfural (kinetic expressions derived from industrial data)

The values of the kinetic parameters are summarized in Table 2.

		1386	4.92
Lactoperoxidase		7110	3.22
Hydroxymethylfurfural	9.98*^16	129684

Bacterial inactivation

Fouling is an important aspect in the operation of heat exchangers for the thermal treatment of milk and milk products, but the primary aim of the process is to obtain a safe product. When trying to  do a model based optimization of the heat treatment configuration, these aspects should also be included in the model to be utilized. According to Kessler (1988) the bacteriological quality of UHT milk can be expressed in terms of destruction of Bacillus stearothermophilus spores.

The parameters can be found in Table 3.

		345400	101.15

Heat exchangers

Heat exchangers are modelled as a 1-D plug systems. For 7-pipes systems it is within reason to assume that the conditions in all pipes are equal so that it is sufficient to simulate the behavior in one pipe and extrapolate the results for the rest. The model equations are constructed by considering the continuity equation for the different components (, , and deposit) over an infinite small volume of the heat exchanger. The continuity equations are then combined with the energy equation over the same volume. The result is a set of coupled partial differential equations (equation 9).

With boundary conditions:

Equations 9 and 10 are set up for a countercurrent type of heat exchanger, as is the case here for the heat exchange subsystems present. For 7-pipes systems, the inlet flow is divided by 7 and is fed to equation 9, and the outlet is then recombined.

Homogenization

The homogenization step is modelled as a small linear increase in temperature. All reactions are included during this small time-temperature interval.

Holders

Holders are modelled as tubular, isothermal systems. Because holders usually have a relative large ratio, it is assumed that no significant fouling takes place. Obviously, the bacteriological and chemical reactions continue to take place.

Integration

Integration of the system partial differential equations is done according to a finite differences scheme with the time dimension as discredits coordinate. For every time step the remaining set of ordinary differential equations is integrated with the NDSolve routine from Mathematica 4, which switches between a non­stiff Adams method and a stiff Gear method.

For every heat exchanger subsystem, all of the countercurrent flow type, the calculation of in- and outlet temperature of the heat exchange medium (water or steam) leads to a 2-point boundary value problem. This is solved iteratively using the shooting method (Press et al., 1992). The heat regeneration couplings in the system lead to comparable problems. For example, the outer-loop coupling has as consequence that the temperature is not known at the first integration of the system. This means that a guess will be needed, after which the system can be integrated. After calculating the system a new value of will be found. If this new value is different from the first guess, the procedure has to be redone. This is in fact done based on an iterative algorithm for root solving (in this case, the secant method is used). The same holds true for the inner loop. The rigorous procedure should thus have to solve the inner loop for every integration of the  outer loop. Taking into account that simulation of one batch period takes ± 500 seconds on a Windows NT workstation, this is in fact not realistic for optimization purposes. There is another problem with the inner heat regeneration loop. The exact configuration is not known, especially the flow rate of the recirculating water and the cooling capacity of the intermediate cooler are important for simulation purposes. It was therefore decided to uncouple the inner heat regeneration loop for now. However, it is planned to include this loop in the near future.

Pressure drop

Fouling causes the hydraulic diameter to decrease during operation. This will cause the pressure drop over the system to rise, up to the point that normal operation is not possible anymore. Obviously, pressure drop is an important parameter and will be considered as the determining factor in comparing different system lay-outs and process conditions. It is calculated using the classic Darcy equation for circular pipes (Perry & Green, 1984) where pressure drop Δp is determined by friction losses and kinetic energy:

The friction factor can be calculated with the Colebrook equation for turbulent flow ( which is definitely the case here):

The surface roughness ε is dependent on the material of the tube. For stainless steel, a value of 0.1 mm can be assumed (Perry & Green, 1984). Several parameters for these equations are temperature or time-dependent. The temperature dependency is taken into account by using an average value for each system module. The time dependency is accounted for in the same manner as for the differential equation integration. A correction factor of 2.5 was used to take deviances from non-ideality into account (couplings, bendings, etc...). With these parameter values, actual pressure drops could be simulated fairly accurate.

Model input

Physico-chemical input

Most of the physico-chemical parameters can be calculated from tabulated data or empirical correlations, presented in Table 4. They were all found in Perry & Green (1984), unless stated otherwise. The value for the density of the fouling layer was calculated from data mentioned in Delplace & Leuliet (1995). is an important parameter because it is partly responsible for the dynamic aspect of the model and it seems probable that its value is not constant as function of process times and location in the heat exchanger. One of the priorities in fouling research should be to establish the time-place dependency of . The details of the calculation as used in this report can be found in the Appendix.

Initial and operating conditions

The operating and initial conditions for the system are summarized in Table 5 and 6. Because of the heat balance problems heat exchanger 2 is basically not used. The initial concentration of βlg and lactulose is in the same order of magnitude as determined for input samples to the actual treatment plant.

	6.11		6.11	0.873	0.873	6.11	0.873
	0.053	0.053	0.085	0.02	0.02	0.057	0.02
	180	120.86	32.14	39.32	38.65	2.39	78.5
	65	43.5	29.8	14.1	13.9	0.99	28.3

	0.0015
	4
	1000
	5
	5
	0.01
	0

Simulation

The system was simulated for a total batch period of 16 hours, which is quite long.

Time-temperature

The simulated time-temperature profile together with the available measured data is shown in Figure 2. It is clear that the model is capable of simulating the temperature quite accurate. The influence of fouling on the temperature is not very pronounced and would not be visible in the Figure. This is due to the fact that the UHT-temperature is imposed, the steam temperature (or steam flow rate) is adapted to the actual product temperature like it is done in practise.

Fouling

The fouling layer thickness as function of the position in the system at the middle (8 h) and end of the production batch (16 h) is shown in Figure 3.  As stated before, it was assumed that no fouling occurs in the homogenization and holder sections. The thickest fouling layer is clearly formed at a process time of about 100 seconds. This is in heat exchanger 3, after the homogenization step. However, this picture is somewhat misleading because heat exchanger 3 is single pipe system with a larger diameter than the 7 pipe systems used for most of the other heat exchanger subsystems. Plotting the ratio of fouling thickness on internal tube diameter gives another view (Figure 4). Looking at this, the fouling production in heat exchangers 4 and 5 are more important than one would conclude based on absolute fouling layer thickness. The relative pressure drop changes (defined as ) for each heat exchanger after the production batch (Table 7) support the view from Figure 4. This Table is important because it summarizes the pressure drop evolution during the production. Optimization strategies will be based on this evolution. The total amount of fouling predicted at the end of the batch in terms of mass deposited is 0.56 kg.

HE1	0.007058428099453005`
HOM	0.011079967272600265`
HE2	1.97872811786605`
HO1	0.0486069284646855`
HE3	2.354605857195773`
HE4	2.261359211673363`
HO2	0.011079905402961361`
HE5	1.2945513177430188`
	1.6814161394989318`

β-lactoglobulin

Figure 5 shows the evolution of the 3 components of βlg as defined in the model of de Jong (see equation 1). As long as there is unfolded β-lactoglobulin present, fouling will occur. Optimization strategies should therefore concentrate on those configurations which favour the transition from unfolded to aggregated β-lactoglobulin which is considered inert with respect to fouling.

Spores destruction

Figure 6 shows the destruction of Bacillus stearothermophilus spores in the system. The evolution at the end of the batch is plotted in dashed lines, but the difference is almost too small to be visible in the figure. The inactivation occurs in a very sharp front, when the temperature gets high enough. However, the system seems to be somewhat overdesigned with respect to spores inactivation. At the beginning of the batch, starting from a (high) initial concentration of 1000 , the outlet concentration reduces to , corresponding to a log reduction of 9.65. This reduction is still 9.47 at the end of the batch. Consequently, the influence of fouling on microbiological safety seems only minimal.

Other components

Hydroxymethylfurfural

The formation of HMF as function of the residence time in the system is shown in Figure 7. The difference between the beginning and end of the production batch is 0.93 which cannot be distinguished in a figure.

Lactulose

Lactulose formation roughly shows the same trend as HMF formation (Figure 8). The difference between beginning and end of production batch is 3.53 .

Alkalic phosphatase

As could be expected, alkalic phosphatase is quite fast destroyed  (Figure 9). The difference between beginning and end concentration is not significant, which could be expected as the reaction takes place in the lower temperature range.

Lactoperoxidase

The evolution of lactoperoxidase shows the same trends as the alkalic phosphatase (Figure 10).

Influence of mass flow

The previous simulations used the maximum possible mass flow for the system: 22000 . In practice the mass flow can vary to values as low as 7000 , so it would be interesting to see what the model output is for other mass flow values. The evolution of the outlet concentration of HMF and lactulose are shown in Figure 11 and 12, respectively. The mass flow has a significant impact on these concentrations, as could be expected. The shape of change is very similar. Obviously, the bacterial destruction is improving with decreasing mass flow. The pressure drop change at is only 0.69 % (comparing with 1.68 % for the basic case, see Table 7). Overall, the pressure drop is lower due to the lower velocity. So, from the viewpoint of bacterial safety and pressure drop increase, lower mass flow is beneficial. It is not beneficial concerning the formation of lactulose and HMF.

Optimization

The system was optimized starting from two different approaches. The first approach was to consider the given system as fixed in the sense that no additional holders can be added. In the second approach an additional holder was added before heat exchanger 3. This is inspired by the fact that reaches a maximum in heat exchanger 3 (Figure 3).

Approach 1

Optimization requires the choice of experimental degrees of freedom. Three were selected: the UHT temperature and the holding time at UHT temperature . can vary continuously between certain boundaries, is limited to 2 values, 1 s and 15 s respectively. Optimization also requires the definition of a cost function. Ideally, this function should comprise elements of process economics, microbiological safety and chemical product quality, each with its own carefully chosen weight factor. The definition of such a cost function is not a trivial task. For the time being, in order not to complicate the calculations too much, only fouling and microbiological safety are included in the cost function. Future research will focus on expanding the cost function definition. The aim is to minimize fouling as much as possible while maintaining a safe product. The concentration of Bacillus stearothermophilus spores at the outlet of the system should therefore be as low as possible. Consequently, a large (exponential) weight is connected to this factor. Fouling is added to this cost function as the maximum pressure drop, as was explained before.

Assuming that we can freely adapt 3 operational parameters , and the residence time in holder 1 , the result clearly point at that fouling reduction can be achieved by as long as possible residence times in both holders and a lower temperature . The corresponding optimal residence times were limited by self-imposed upper limits ( and ). This confirms that it seems more beneficial to operate at a lower temperature at a longer time, in terms of fouling. The bacterial inactivation is still guaranteed at these values, as is evident in Table 8.

		136.85	133.4
	1	15	42
	30	30	92.9
	1.68	1.41	0.85
	9.65	8.54	8.10

Approach 2

An additional holder before heat exchanger 3 was added to the model. Additionally, the inlet temperature for this additional holder was assumed to be adjustable between certain limits. In practice, this would require heating, e.g. in heat exchanger 2. Figure 13 shows the evolution of the pressure drop change as function of additional holding time and inlet temperature. Obviously, the temperature in the additional holder has a more pronounced influence than the residence time . In fact, the optimization roughly coincides with a longer residence time in holder 1. The influence of the residence time in holder 1 is shown in Figure 14, where it can be seen that the reduction can be quite pronounced. Thus, instead of inserting an additional holder, there can be quite some benefit by simply prolonging . Although Figure 13 shows that there is an optimal temperature for an additional holder (± 90 °C), the margins are very small.

Combining approach 1 and 2 () leads to a fouling reduction of 51 %.

Influence on the formation of other components

Adapting process conditions can be expected to have an influence on the formation of some components. Table 9 show the predicted formation of HMF and lactulose for the different simulated setups. Obviously, the UHT-conditions have a larger influence than the additional holder, despite the long residence time.

Basic		329.56
		426.99
	138.42	615.37
		338.10
	98.98	436.17

Conclusion

The model simulations show that it is possible to reduce fouling, as expressed by pressure drop change, using only the available degrees of freedom by setting the UHT temperature lower and the residence time higher. Extending the intermediate holder 1 has also a positive effect. The latter also minimizes the extra formation of lactulose and HMF, whereas a longer UHT residence time has a profound negative effect. On the other hand, lowering the reduces the energy consumption of the process. Varying mass flow has also a strong impact on the formation of these components. Due to the lower velocity, the pressure drop change is less with lower mass flow, and bacterial inactivation is higher.

Symbols

Plain text

C	Concentration
C
D	diameter	m
E		kW
F
	Fouling
k
M
N
p	Pressure	Pa
R
r
T	Temperature	K
t	Time	s
U
v	Velocity
z		m

Subscripts

AHO
AP
bst
CO	Cooler
HE
HMF	Hydroxymethylfurfural
HO	Holder
lac	Lactulose
LPO	Lactoperoxidase
wall

Greek

δ	Thickness	m
Δ	Difference
λ
ρ	Density
τ		s
φ

References

Delplace F. & Leuliet J.C. (1995). Modelling fouling of a plate heat exchanger with different flow arrangements by whey protein solutions. Trans. IChemE Part C. 73, 112-120.
de Jong P. (1996). Modelling and optimization of thermal processes in the dairy industry. NIZO Research Report V341, Ede, The Netherlands.
Froment G.F. & Bisschoff K.B. (1990). Chemical reactor analysis and design. Wiley, New York, USA.
Kessler H.G. (1988). Lebensmittel-und Bioverfahrenstechnik; Molkereitechnologie. Kessler Verlag, Weihenstephan, Germany.
Peri C., Pagliarini E.& Pierucci S. (1988). A study on optimizing heat treatment of milk. I. Pasteurization. Milchwissenschaft, 43, 636-639.
Perry R.H. & Green D. (1984). Perry's chemical engineers' handbook. McGraw-Hill International Editions, New York, USA.
Press W.H., Teukolsky S.A., Vetterling W.T. & Flannery B.P. (1992). Numerical recipes in C. The art of scientific computing. 2 Edition. Cambridge University Press. Cambridge, UK.
Walstra P., Jeurink T.J., Noomen A., Jellema A. & Van Boekel M.A.J.S. (1999). Dairy Technologies. Principles of milk properties and processing. Marcel Dekker, New York, USA.

Appendix 1: Calculation of the fouling layer density

The fouling layer density was estimated based on data presented in a paper from Delplace & Leuliet (1995). They measured the apparent thermal conductivity as function of process time for a certain heat exchanger configuration. They deduced the following relationship between and :

All the necessary data needed to calculate is present in the paper (Tables and Figures). Table 10 lists the values used.

N		6
t

Substituting these values into equation 13 leads to a value of 940 for . This is the value for the density of the fouling material, without taking the porosity into account. According to Walstra et al. (1999), the porosity in the fouling layer is in the order of magnitude of 55 %. This means that the apparent density of the fouling layer, , would be 380 .
One has to keep in mind that these values must be considered approximate. The problem with (or ) is that its value has a profound impact on the simulation results. Figure 15 shows the fouling mass deposited as function of , where the arrow indicates the value used for the calculations in this report. The importance of on the model output is obvious. Unfortunately, the real value is very difficult to determine. As long as there is significant uncertainty on this value, the results of any simulation exercise should be interpreted with caution, and should be looked at as pointing towards certain trends whereas the absolute numbers should be given less importance.

Appendix 2: Calculation of energy consumption

The energy used to heat the milk will heavily influence the economy of any heat treatment process. Therefore, an approximate calculation to assess the energy consumption of the system was worked out. Approximative, because the inner heat recuperation cycle is not included in the simulation. Referring to Figure 1, the modules where external heat is added to the system are heating 2 and the UHT heat exchanger. Heat exchanger 2 is basically not used, so only heat exchanger 5 remains. At every time step, the amount of energy added to the system per unit of length follows directly from the energy equation 9. It has to be remarked that the cooling costs for the inner -cycle are not included. The overall energy added to the system is calculated by integrating the energy consumption per unit of length over the total length of the heat exchanger, which is then averaged over the process batch time. The results are summarized in Table 11. Obviously, lowering the UHT temperature has a beneficial effect on the energy consumed.

Basic	57.0
	31.3
	24.5
	57.3
	32.2