- Email: [email protected]
Prediction of freezing and thawing times for multi-dimensional shapes by numerical methods
Prediction of freezing and thawing times for multi-dimensional shapes by numerical methods D. J. Cleland, A. C. Cleland, R. L. Earle and S. J. Byrne* Massey University, Palmerston North, New Zealand Received 21 April 1986
An assessment of the accuracy of numerical methods used in the prediction of freezing and thawmg times wa~ made using a comprehensive set of freezing and thawing data for both regular and irregular multidimensional shapes. For regular shapes, a finite difference method gave accurate predictions with reasonable computation costs. Predictions for two finite element method formulations were not always accurate. This was due to practical constraints on the computation costs which meant that time and spatial grids could not always be made sufficiently fine to ensure that the prediction method uncertainty was insignificant compared with the other sources of imprecision. Guidelines are suggested for the use of the finite element method as a freezing or thawing time predictor. These should ensure that the prediction method error is small while keeping the computation costs reasonable.
(Keywords:freezing;thawing;computerprograms)
PrOvision par des m6thodes num riques des temps de congelation et de dbcong61ation pour des formes multidimensionnelles Une Ovaluation de la prbcision de mbthodes numOriques utilisOes pour la provtslon des temps de con qe.lation et de dOeonoblation a btb effeetuOe en utilisant un ensemble complet de donnOes sur la conoblation et la dOcongblatton de formes multidimensionnelles rOguliOres et irrOguliOres. Pour les formes rbguliOres la mOthode des difference~s finies a donnO des provisions prbcises avec des cofits de calcul raisonnables. Les provisions pour deux formules de la mOthode des Olbments finis n'Otaient pas toujours prOcises. C'Otait dfi h des contraintes pratiques sur les cofit~ de calcul, ce qui signifiait que les 9rilles de temps et d'espace ne pouvaient pas toujours ~tre suffisamment fine,s pour s'assurer que rincertitude de la mOthode de provision Otait peu importante en comparaison avec d'autres sources d'impr~cisions. On propose des directives h utiliser dans la mOthode des OlOments finis comme moyen de provision du temps de conoOlation ou de dOcongOlat ion. Avec ces directives l'erreur de la mOthode de provision est faible tout en maintenant des cofits de calcul raisonnables.
(Mots cles: cong~lation; d~cong~latlon; programmes d'ordinateur
Nomenclature Bi D 2D, 3D FDM FEM h k SD
Biot number = hD/k Object thickness or characteristic dimension (m) Two-dimensional, three-dimensional Finite difference method Finite element method Surface heat transfer coefficient ( W m - 2 oC -1 ) Thermal conductivity ( W m -2 '-C 1) Standard deviation of percentage differences
During freezing or thawing of food materials latent heat is released or absorbed over a range of temperatures below O°C. The most physically realistic model for freezing and thawing of foodstuffs is that of transient non-linear heat conduction with temperature-variable thermal properties in the solid object, subject to convective heating or cooling at the object surface1. Without making grossly
* Present address: Auckland University, Auckland, New Zealand 0140-7007/87/010032~)8503.00 Q~ 1987 Butterworth & Co {Publishers} Ltd and IIR
32
Int. d. Refrig. 1987 Vol 10 January
T
Temperature t Ct
Subscripts
a f in x, v, z
Ambient freezing or thawing medium Full unmodified formulation Initial Simplified formulauon In the x, y, or z direction, respectively
simplifying assumptions, the non-linear pamal differential equations representing this model cannot be solved exactly by analytical methods to predict freezing or thawing times. A wide variety of approximate methods have been suggested as solutions 2. These methods fall into two broad groups: numerical methods; and simple formulae. Both groups have advantages and disadvantages so it is expected that food engineers will continue to make use of both for some time to come 1'3 The use of simple formulae is discussed elsewhere 4"5.
Prediction of freezing and thawing times: D. J. Cleland et al. Numerical methods The numerical methods make discrete mathematical approximations to the time and spatial variations defined by the governing partial differential equations for heat conduction. If they are formulated and implemented correctly to reduce numerical truncation and rounding errors, the numerical methods are generally considered the most accurate, reliable and versatile freezing and thawing time prediction methods 1'°. The numerical methods are further subdivided into two major types: finite difference methods; and finite element methods. Both have been extensively used for prediction of freezing or thawing rates and times 2'7'8. Detailed discussion of the advantages and disadvantages of the two methods can be found in previous papers 8"9. The finite difference method (FDM) is limited by practical considerations to use for one-dimensional and regular multi-dimensional shapes where evenly spaced grids can he used, while the finite element method (FEM) can more easily take account of irregular multi-dimensional geometries. Although the spatial discretization is different for the FEM it is usual to employ a FDM scheme for the time progression in transient problems. Bonacina and Comini l° and Cleland and Earle 1 considered different F D M schemes for phase change problems and concluded that the implicit, three time level, Lees scheme 11 was most accurate. Formulations of both the F D M and FEM based on the Lees scheme have been developed and implemented to predict freezing and thawing for a range of geometries s'9' 12.13,14 It is well known that with a very fast computer, unlimited computer memory and appropriate software both time and space intervals can be made sufficiently small that the F D M and FEM will give accurate predictions for multi-dimensional shapes, whether they are irregular or not. But published results indicate that few researchers, let alone industrial food engineers, have such computational resources routinely available. Therefore, it cannot be assumed that the F D M and FEM are universally accurate predictors of freezing and thawing times in practical situations. The size of the available computer memory limits the refinement of the spatial grid that can be used. The sizes of both the time and space intervals are restricted by the computer processing speed and cost. In some situations these limits can be overcome. For example, a F D M used with an 8 × 8 × 8 evenly spaced node grid that required less than 0.5 megabytes of memory and took less than 1000s process time on a Prime 750 computer gave good predictions for freezing of rectangular brick shapes. By comparison, a FEM using a 6 × 6 × 6 evenly spaced nodal grid required about 2 megabytes of memory and 30000s of process time, and yet gave barely adequate prediction accuracy for the same problem 9. Because of such difficulties, modifications to the basic F D M and/or FEM formulations have been developed either to allow the use of a larger time step while maintaining the same prediction accuracy or to reduce the computation effort per time step 9'12't5"1°'17"18'19. With all of these modifications there is a trade-off between the loss or gain in accuracy and the savings or increase in computation costs. In practice, uncertainties in thermal property and boundary conditions data are often substantial, so that provided numerical approximation errors are kept small compared with the physical data
uncertainties, no significant change in prediction accuracy will result. For multi-dimensional problems, users may resort to using coarse spatial grids to save computation time. This can lead to prediction inaccuracy as the spatial approximations may become limiting in accuracy, even when one of the modified methods is used. Another major disadvantage of numerical methods compared with simple formulae is that they require detailed thermal property data over the whole temperature range. There are a number of methods (see, for example, Heldman 2°) to predict thermal properties over a wide range of temperatures from only a few basic measurements, such as moisture content and the initial freezing temperature, so that reasonably accurate prediction of thermal property data is usually possible. Calculation of a heat balance 1'9 as part of the numerical method implementation is one of the best ways to check the accuracy of all the numerical approximations made by the numerical methods 1. The calculated heat flow through the surface of the object should equal the change in the object's internal enthalpy. Close agreement (within ~ 2 °o) should be achieved for a numerical method to be considered both accurate and valid ". Assessment of prediction accuracy o/ numerical methods For problems that can be considered in one dimension (slabs, infinite cylinders and spheres) both the F D M and the FEM are easily implemented, require minimal computational effort and give comparable prediction accuracy. By comparison with experimental data for freezing 3"9'2x and thawing 5 of typical foods, three time level versions of both methods have been proved to be sufficiently accurate predictors of both freezing and thawing times. For freezing of two- and three-dimensional regular shapes (infinite rectangular rods, finite cylinders and rectangular bricks) the accuracy of the Lees scheme F D M has been verified4' 13.22.23. As expected for regular shapes, a full, three time level finite element method formulation (FEMf) was found to give comparable results to the F D M if sufficiently fine time and space grids were used 9. A simplified finite element method formulation (FEMs) 9 was found to be nearly as accurate as the FEMf formulation and gave large savings in computation time. However, it was not clear how the FEM, should be used as, in all cases tested, the spatial grids were reasonably fine and the time steps were controlled so that good heat balances were maintained. For phase change in multi-dimensional irregular shaped objects, comparisons of FEM predictions with experimental data have been made for freezing of roads24, thawing of fish 25, freezing of ice 26'27 and freezing of beefz8'29, all in two dimensions, and for freezing offish 3° in three dimensions. These comparisons were all limited in scope, as they were made for only one or two experimental runs each. Research aims In summary, the accuracy of the three time level F D M has been proved for regular shapes within reasonable computational cost limitations. Although it was expected that the FEM would give accurate prediction of phase change this testing has not been comprehensive,
Rev. Int. Froid 1987 Vol 10 Janvier
33
Prediction of freezing and thawing times: D. J. Cleland et al. Table 1 Percentage differences between the experimental freezing and thawing Umes for Tylose mulU-dlmensumal ,hapes and the freezing and thawing times calculated by numerical methods Tableau 1 Difference en pourcent entre les temps de congelation et de decon.qi'latlon expbrimentaux pour des h~rme~ muhuhmensumnelle~ de l-vlow eL te~ temps de congblation et de dbcongblation calculbs par des mbthodes numbrtque~ Mean (%)
SD (",,1
MmJmum {',,)
Maxtmum ~" ~
Shape
Data 23, 31
Method
Rectangular bricks
Freezing (72 runs)
FDM FEM,
4,4 L9
6.6 7 1
2u 24
Thawing (68 runs)
FDM FEM s
0.6 b3
3 1 34
(, ,~ 12
'~
Freezing and thawing (83 runs)
FEM~ FEM,
0.3 42
5 I 6.1
I! ~ ')~
~l : 2) ~'
Freezmg (42 runs) Thawing (41 runs)
FEMf FEM,
10 26
45 4.8
II ~ ,1
'~, I ? ~)
FEMr FEMs
15 58
55 70
(,
i~ ?1 ,)
Two-&mensmnal irregular shapes
especially for multi-dimensional irregular shapes, and where there are computation cost restrictions. Therefore the aims of this paper are: 1. to assess a three time level finite element method (FEM09 for the prediction of both freezing and thawing times for multi-dimensional irregular shapes by testing against experimental data subject to practical computational cost constraints; 2. to assess the value of the simplified finite element method formulation (FEM,) as proposed by Cleland et al. 9 as an alternative to save computation costs; 3. to formulate guidelines for the use of the FEM for phase change problems that would help to ensure accurate predictions where there are computation cost limitations; and 4. to use data for thawing of rectangular brick shapes to verify further the accuracy of the numerical methods for regular multi-dimensional shapes. Results and discussion
A large set of experimental data was collected for this assessment of the numerical methods. Details of experimental techniques and the experimental results are reported by Cleland et al.31. Assessment of the prediction accuracy of the numerical methods was made by direct comparison of the numerically calculated times for the thermodynamic centre temperature of the objects to reach - 1 0 ° C for freezing, or 0°C for thawing, with the experimentally measured times. The percentage differences were calculated from: 3/0 Difference =
Predicted time-Experimental time 100 ×--Experimental time 1 (1)
In all cases the thermal property data used for the numerical method predictions were those given by Cleland et al. s and at least 1000 time steps were used in each numerical calculation. Prediction of freezing and thawing times for rectangular bricks Table I summarizes the prediction accuracy of the
34
Int. J. Refrig. 1987 Vol 10 January
,~ s
-, ,
F D M 23 and the simphfied finite element method formulation, F E M f , compared with experimental data for freezing23 and thawing 31 of rectangular brick shapes of Tylose. The full finite element method formulation. FEMf 9, was not used for the full data set because computation costs were large with the grids used. If the FEMf had been used, the results would be expected to match the F D M results 9 closely. Because of symmetry only an octant of each rectangular brick was modelled. Within this octant an evenly spaced 9 x 9 x 9 node grid was used in the F D M calculations. The FEM calculations used a 7 ~ 7 x 7 node grid with 216 evenly sized, eight-node, linear isoparametric, brickshaped elements. A finer grid was not practicable due to limitations on computation cost and computer memory size. For all except a few runs the calculated heat balance/'9 agreed to within 2 o . Computer process times on a Prime 750 computer were typically of the order of 1000s per run for the F D M and 10000s per run for the FEM,. The FEMr was about five times slower than the FEM,. The freezing results have been discussed previously ". The high spread of percentage differences and the overall under-prediction was attributed to freezing rate affecting the precision of the thermal property data rather than to uncertainty in the experimental data or imprecision due to the numerical methods 13. Thawing is inherently slower than freezing, so the same effect was not expected for the thawing data. For thawing, Table 1 shows that the FEMs formulation had only a slightly greater offset of the mean prediction error than the predictions for the FDM. As was the case for thawing of the simple slab, infinite cylinder and sphere shapes 5, this difference was thought to be due to the relatively crude way in which variations in the thermal properties are incorporated into the FEMs. The relatively coarse nodal grid used for the three-dimensional geometry accentuated this problem, so the difference between the F D M and FEM, is greater for rectangular bricks than for the simple one-dimensional shapes. Predicted temperature profiles agreed closely with experimentally measured profiles, as illustrated in Figure 1. Uncertainty in thermocouple placement was probably a major contributor to the differences between them. Both the F D M and FEM, formulations predicted oscillating
Prediction of freezing and thawing times. D. J. Cleland et al. 40 I
I E
0
-20 0
I
I
I
I
1
2
5
4
Time (hrs) Figure I Typical temperature l,ersu.~ time profiles for thawing a rectangular brick of Tylose. Run 7 (Table 1 of Reference 31) D~= 0.104m; D,.=0.125m: D~=0.152m; h = 4 1 . 0 W m -2 'C - i , Ta= 45 0 C; T~. = - 12.0 C. , Experimentally measured: . predicted by F D M : . predicted by F E M s. 1. Temperature at the centre of the :,,z surface: 2. temperature at the geometric centre Figure 1 Dtagramme type tempbrature temps pour la dbcong~latton d'un bloc parallel epipbdtque de 7].lose. Essat 7 (le tableau I de la r~ference 31). Dx=O,lO4m. D y - O , 1 2 5 m . Dz=O,152m; h = 4 1 . O W m 2C 1. T , = 4 5 . 0 C. T,, = - 12.0 C --. Me~ur~ e,:p~rimentalement. . prbt u par F D M. . prbt,u par F E M , 1. Tempbrature au centre de la ~urface x.z. 2, temperature au centre geomi, trique
surface temperatures due to the discrete nature of the numerical methods. For the FEM~ formulation this was magnified by the coarse grid used for the predictions 4. Use of the three time level enthalpy method ~8 or the other modifications to the method to improve estimation of the thermal properties in the phase change temperature range ~2'15 ~ may help reduce this type of oscillation. The differences between predicted and measured thawing times were not statistically correlated with any independent variable in the experiments. Further, the 95 °'o confidence bounds for the thawing time predictions of the numerical methods are similar in magnitude to the estimated experimental uncertainty bounds of + 5.4 °,o31. This suggests that for the F D M and FEMr formulation, most of the difference between the predictions and experimental values was due to experimental uncertainty and not imprecision resulting from the implementation of the numerical methods themselves. The only significant inaccuracy arising from the numerical methods was that introduced by the crude thermal property estimation in the FEM, formulation, thus indicating that care should be exercised in using the simplified formulation with coarse grids.
Predictions fi~r two-dimensional irregular shapes The experimental freezing and thawing times for twodimensional irregular Tylose shapes are given in Table 2 of Reference 31. The grids shown in Figures 1-7 of Reference 31 were used to make predictions by both FEM formulations. Eight-node, quadratic, isoparametric, curved rectangular elements were used. The grids prepared ensured that at least five elements (11 nodes) existed between the thermodynamic centre and the outside surface of each object. The FEMf formulation
gave satisfactory heat balances for all the runs, but the FEM s formulation did not always achieve heat balances within the 2 ° o error bounds, although computation costs were reduced by about 80°4, for the same time step size. A typical calculation for the FEMf required about 2 megabytes of memory and took about 50000s process time on a Prime 750 computer. The results, summarized in Table 1, show that both FEM formulations gave accurate predictions of freezing and thawing times (within the tolerances allowed for experimental error) for all the shapes tested. The 95 ° o confidence bounds for the percentage differences of the FEMf formulation of - 9 . 7 to 10.3 ° o are similar to the estimated experimental error of + 8.9 ",,31 indicating that most of the difference was explainable by the experimental uncertainty. Again, the FEM, formulation did introduce some extra inaccuracy due to the cruder thermal property estimation, as indicated by the slightly higher spread and mean offset of the percentage differences. The freezing data were under-predicted compared with the thawing data by both formulations. This was attributed to the freezing rate effect t3 rather than numerical approximation problems because good heat balances were consistently achieved for the FEMf formulation. Typical predicted temperature profiles (Figures 2 and 3) compared well with those measured experimentally. Again, differences between them are probably due in part to uncertainty in thermocouple placement. Also, the possibility of undesired variations m the local surface heat transfer coefficient for these experiments 4'5 meant that measured temperatures at points on, or near, the surface had significant uncertainties associated with them. The surface temperatures predicted by the FEMs formulation
,0 7 20-
oQ_
E -20 i
- 4 0 ~ - -
0
- - F
2
-
-
T . . . . . .
I
4
6
--
---q
8
Time (hrs) Figure 2
Typical temperature rersu~ time profiles for freezing a Tylose two-dimensional irregular shape. Run 29 (Table 2 of Reference 31), s h a p e n o 3. D = 0 . 1 3 7 m ; h = 2 8 . 0 W m - 2 C 1; T a = _ 3 0 . 5 o C , Tm= 32.8 C. - - , Experimentally measured: , predicted by FEMf, -- , predicted by F E M s. The object shape and thermocouple positions 2, 3 and 5 are shown in Figure 3 of Reference 31 Figure 2 Dlagramme type tempbrature temp,s pour la congblation d'une forme irr#guh#re de Tylose gt deux dimensions. Essal 29 (tableau 2 de rbference 31). Forme no. 3 D = O , 1 3 7 m . h = 2 8 , 0 W m 2 'C I T, = - 30.5'~C. T~ = 32,8 'C - - , Mesur~ expbrimentalement. -. prbvu par F E M f: -. prei,u par F E M s. La forme de l'objet et les emplacements du thermocouple sont indtqukes .figure 3 de la rbference 31
Rev. Int. Froid 1987 Vol 10 Janvier
35
Prediction of freezing and thawing times." D. J. Cleland et al. heat balances were not achieved for all of the predictions. mainly because spatial grids were coarse for all of the calculations with three-dimensional grids. The size of the available computer memory, the cost of the computer processing and/or the large data preparation time meant than finer grids could not be used. The grid used for the pyramid shape is shown m Figure 8 of Reference 31. It is basically the same 7 x 7 × 7 node grid used for the rectangular brick predictions, but is distorted to the pyramid dimensions. Because of the three-dimensional irregularity a full quadrant of the shape was modelled. Consequently the grid was coarse. with only three linear elements between surface and centre along the height of the pyramid. Computation costs were the same as for the rectangular brick predictions. The predictions have similar spread to those for freezing and thawing of rectangular bricks but the mean was offset by about 10°o . The low standard deviation of the prediction errors suggests that the source of this error is systematic, not random. Two possible reasons are: the coarse grid leads to significant prediction inaccuracy: and the assumption that h was the same for the pyramid and the reactangular brick shapes 3~ may not hold. Both reasons probably contributed. In the numerical calculations, freezing and thawing were considered complete when the slowest cooling or heating node reached the desired final temperature. If the true thermodynamic centre of the object does not exactly coincide with the slowest node position at this time, then the object will not be completely frozen or thawed. Under-prediction of the time for the complete process will result from this effect. The effect is smallest for fine finite element grids, but for coarse, poorly positioned grids the extra time could be substantial. Further, linear elements magnify the problem because they use linear approximation functions between nodes. The true temperature profile and any local maximum or minimum values of the temperature within each element may not be properly modelled. Quadratic and higher-order elements allow more accurate approximations to the temperature profile to be made and therefore often give more accurate
oscillated significantly compared with the predictions for the FEMf formulation. Again, this difference was attributed to the crude thermal property estimation used in the simple formulation, which was accentuated by the relatively small number of quadratic elements in the grids used for the predictions.
Predictions for three-dimensional irregular shapes Table 2 summarizes the percentage difference between experimental freezing and thawing times for some threedimensional irregular shapes (Table 3 of Reference 31) and the predictions by the two FEM formulations. Good
10
? v t-
Ill
t:l.
E
-10-
-20 --~ I
-30-~--0
I
[. . . . .
2
4
I
I
I
6
8
10
.....
7
12
Time (hrs) Figure 3 Typical t e m p e r a t u r e versus t i m e profiles for t h a w i n g a Tylose t w o - d i m e n s i o n a l irregular shape. Run 31 (Table 2 of Reference 31 ), s h a p e no. 4: D = 0 . 0 9 9 m : h = 2 8 . 0 W m - 2 " C - ~ : T a = 2 0 . 1 ° C : Tin= - 21.5'-C. - - , Experimentally measured: , predicted by F E M f , -, predicted by F E M s. The object shape a n d t h e r m o c o u p l e positions 2, 3 a n d 6 are s h o w n in F i g u r e 4 of Reference 31 F i g u r e 3 Dlagramme type tempbrature temps pour la dbconoblatum
d'une forme irrboulibre de Tylose h 2 dimensions Essat 31 (tableau 2 de rbference 31), forme no. 4 D=O,O99m. h = 2 8 , 0 W m -2 C -1 T = 2 0 , 1 ~ C . T~ = - - ~1,,~ C , Mesurb experirnentulement, prbvu par F E M y. -- -, prbvu par FEM,. La forme de l'objet et les emplacements du thermocouple sont indlqubs figure 4 de rbfferenee 31
Table 2 Percentage differences between the e x p e r i m e n t a l freezing and t h a w i n g u m e s for l y l o s e t h r e e - d i m e n s i o n a l irregular shapes and the freezing a n d t h a w i n g times c a l c u l a t e d by n u m e r i c a l m e t h o d s T a b l e a u 2 Difference en pourcent entre les temp~ de congblatlon et de dbcon,qielatton expermwntaux pour des ]brines muh Mlmen~ionnelles lrrequhere~s dc
Tylose et les temps de congblation et de dbcongblation calculbs par des mbthodes numbriques
D a t a 31
Grid
Method
Mean I ",,)
Pyramid {6 runsl
3D
FEMf FEM,
-- 10.3 - 14 5
6.7 4 8
IS 2 _t1,_ ~ ~
- 20 8q
3D
FEMf FEM~
- 5.5 --I 2
6 8 12 7
- t2 I lt~ I
5.3 13 9
2D
FEMf FEM s
0.6 7.7
70 8.3
,~ 7 J6
13.9 18 2
ID
FEMf FEM s
0.4 2,6
67 65
,' v 4.9
9 6 11 6
3D
FEMf FEM s
5.2 9 8
13 0 14 9
-- 12.9 - ~) 4
25.1 28.8
2D
FEMf FEMs
3,7 16.4
13.8 18 7
- lh 1 4 1
24 7 48 4
3D
FEMf FEM s
35.3 33.5
89 16.3
20 8 8.6
44 2 58.4
Sphere (6 runs)
Egg
(6 runs)
Fish (6 runs)
36
Int. J. Refrig. 1987 Vol 10 January
SD I",)
Minimum 1" ~
Maximum I ",,t
Prediction of freezing and thawing times. D. J. Cleland et al. predictions. FEM calculations for the pyramid shape were affected to some extent in the manner just described as a reasonably coarse grid of linear elements was used, The estimate of the position of the thermodynamic centre was not accurate so the 'centre' node position was not very close to the actual thermodynamic centre. The amount of offset was greater for the freezing runs where the Biot number, Bi, was low, suggesting, according to the criteria of Cleland and Earle ~, that the second reason (systematic error in h) was also an important contributor. The other three shapes were a sphere, an egg and a fish shape. The sphere is a regular shape and can be modelled by one- or two-dimensional grids rather than threedimensionally if axisymmetrical formulations of the FE M are used. Three different grids were used (Figures 1 and 11 of Reference 31 i: 1. an 11 node, 10 linear element, one-dimensional grid (typical computer process time of 100 s per run on a Prime 750 computer): 2. a two-dimensional grid modelling a quadrant of the cross section (a circle) with 21, eight-node, quadratic, isoparametric elements (typical process time of 5000 s): and 3. a 208 node, three-dimensional grid modelling an octant of the sphere with 27 quadratic, isoparametric elements (typical process time of 70 000 s), The three-dimensional grid was coarse, with only three quadratic elements between the surface and the centre. These were arranged in three expanding shells so that the elements were all approximately the same size and none was too distorted from the basic rectangular brick shape. The sphere shape predictions best illustrate the effect of using a coarse grid and the difference between the FEMf and FEM~ formulations. The results generated using both the finer one- and two-dimensional grids with the FEMf formulation were very similar. The results generated using the coarser three-dimensional grid were significantly less accurate. Incorrect implementation of the three-dimensional part of the program was unlikely to be the cause as the program had been extensively tested 9. The position of the thermodynamic centre was known exactly for the sphere, and quadratic elements were used. Therefore, the problem with the coarse grid was less acute than for the pyramid shape predictions and mean underprediction for the coarser three-dimensional grid, compared with the finer one- and two-dimensional grids, was less than the mean under-prediction for the pyramid shape. For the sphere, the FEMs formulation gave comparable results with the formulation, but where the element size was larger in the two- and threedimensional grids the results had a larger variation due to the less representative incorporation of thermal properties into the finite element scheme. For both finite element methods the freezing and thawing time prediction accuracy as well as the prediction of temperature profiles was similar to that achieved for a large set of sphere thawing data 5. The experimental techniques were similar, except for use of the plastic moulding for the shapes, so similar prediction accuracy was expected. The egg shape has one axis of rotational symmetry so it could be modelled by either a two- or three-dimensional grid (Figure 9 of Reference 31). Both the grids were
double versions of the sphere grids, with one half elongated to fit the ovoid shape. The fish shape, although irregular in all three dimensions, was closely approximated by a shape with two planes of symmetry, so only a quadrant was modelled. The grid is shown in Figures l0 and 12 of Reference 31. It is the same basic three-dimensional grid used for the egg shape, but the nodal coordinates, especially for the surface nodes, were modified to fit the actual shape. A typical calculation using the threedimensional egg grid or for the fish shape took 150000s process time on a Prime 750 computer and required 4 megabytes of memory. For the egg and fish shapes it was assumed that the surface heat transfer coefficient was the same as that measured on the moulded plastic spheres 31. The additional inaccuracy introduced by this assumption could not be easily assessed. This extra experimental uncertainty, plus the uncertainty due to difficulties in measuring and describing these shapes mathematically, as well as the approximation of the shape with a coarse FEM grid, meant that poorer agreement was expected for the egg and fish shapes than for the sphere and pyramid shapes. For the egg shape the mean prediction error was close to zero, so the assumption that the surface heat transfer coefficient was the same as it was for the sphere seemed reasonable. The higher standard deviation of the predictions for the egg shape compared with the pyramid and sphere shapes was considered to be mainly the result of the additional uncertainty in the experiments rather than uncertainty introduced in the implementation of the FEM. Ifa finer grid had been used, the FEMf would not be expected to predict the experimental data for the egg much better. The accuracy of temperature profiles for the egg shape (Figure 4), calculated using the three-dimensional grid reflected the accuracy of the freezing and thawing time predictions. Temperature profiles for some surface positions may have been affected by localized variation in h 4. The predicted temperature profiles were affected by the coarse three-dimensional grids used, and therefore displayed some oscillatory behaviour. As expected, this was particularly obvious for the FEM, formulation predictions, for positions near the object surface and at temperatures close to the phase change temperature range 9. For the fish shape the predictions of freezing and thawing times were consistently high by 20-40 °,o but the standard deviation of prediction errors was comparable with that for the egg shape. The offset was probably caused by systematic error in the surface heat transfer coefficient, which may not have been the same as that for the sphere or egg due to: thinner plastic thickness arising during moulding; and/or better thermal contact between the Tylose and the plastic mould for the fish shape than for the sphere and egg shapes. As discussed previously 31, these factors could not be quantitatively assessed as it was not possible to measure h directly for the fish shape. However, measurement of the plastic thickness showed variations of up to 15 O/o of the average total thickness over the surface of the fish shape and differences of up to 10 °/o in the thickness between the fish and sphere shapes. Further, predicted temperature profiles for the fish shape 4 showed trends consistent with a systematic underprediction of h.
Rev. Int. Froid 1987 Vol 10 Janvier
37
Prediction of freezing and thawing times. D. J. Cleland et al. and the two-dimensional irregular shapes. This was attributed to systematic error in the estimate of the surface heat transfer coefficient. The assumption that h was the same for the meat experiments as for the Tylose experiments may not have been valid 3~. Tylose and minced lean beef have different textures, so the contact resistance between the object mould walls and the materials may not have been equal. In view of the small data set, firm conclusions cannot be made about prediction accuracy for the meat data, Howe~er, the general agreement should be noted
Prediction of freezing and thawing data for minced lean beef Freezing and thawing data for minced lean beef rectangular bricks and two-dimensional irregular shapes (Table 4 of Reference 31) were also predicted using the FDM and FEMf. The same grids employed previously were used for the calculations. Table 3 summarizes the results. It was expected that the numerical methods would introduce negligible additional prediction inaccuracy when applied to freezing and thawing of actual foods. However, an increase in thermal property data and experimental uncertainty relative to the Tylose data was expected4. The larger standard deviation of the prediction errors for the bricks compared with the equivalent Tylose results may reflect the difficulties in experimenting with actual food materials. The mean percentage differences for the meat data were higher than for Tylose by about 10% for both the bricks
Guidelines ,[br use of the finite element method Some guidelines for application of the FEM were established by Cleland et al. 9. In the light of the present assessment against experimental data these can be refined. Adherence to these guidelines should lead to predictions in which the error arising from the numerical approximations is negligible compared with experimental and data uncertainty while keeping computation costs as low as possible. The guidelines are'
30-
1. to use sufficiently small nine steps that the heat balance agrees to within 2".~: 2. to use at least 8-10 nodes and at least five elements between the thermodynamic centre and the surface in the FEM grids; 3. to position enough nodes near the estimated position of the thermodynamic centre to ensure that one will ultimately be reasonably close; and 4. to use the full FEM formulation if computational cost limitations allow.
20-.
o
~
-10 -
Of these guidelines, the first two are the most important. In almost all situations, using more refined grids and smaller time steps than those suggested will only give small increases in accuracy and will probably not justify the increase in computation costs
-20 -
-300
1
I
1
I
1
1
2
3
4
5
I. . . .
q
6
7
Conclusions
Time (hrs)
Practical constraints on computation power mean that users of numerical methods, especially the FEM, will need to choose the grid size carefully. The finite difference method used for regular two- and three-dimensional shapes predicted thawing times a~ accurately as it predicted freezing times. Even with computation cost constraints the method could be implemented with negligible inaccuracy arising from the numerical approximations used. The implementations of the FEM used for both regular and irregular multi-dimensional shapes did not lead to significant inaccuracy in predictions provided the
Figure 4 Typical temperature ver,~us time profiles for freezing or thawing a Tylose three-dimensional irregular shape. Run 16 [Table 3 of Reference 31), egg shape: D = 0 A 7 0 m ; h = 5 1 . 4 W m - 2 C ~, T a = - 2 2 . 6 ° C ; Tm=20.6°C. - - , Experimentally measured. . predicted by FEMf; - - - - - , predicted by F E M s. The object shape and thermocouple positions 2, 4 and 6 are shown in Figure 4 of Reference 31. calculated using the three-dimensional grid Figure 4 Diagramme type tempbrature-temps pour la congblatton et la dbcong~lation d'une forme irr~gulibre de Tylose it 3 dimensions. Essai 16 (tableau 3 de rbference 31). Forme de l'oeuf D=O,17Om. h = 5 1 , 4 W m - 2 o c 7, T = _ 2 2 , 6 o c . T =20,6~C. , Me,surb expbrimentalement, - - , prbvu par F E M y, , prbvu par F E M s La forme de l' objet et les emplacements du thermoeouple sont indiqubs figure 4 de la rbference 31, calculbes en utilisant la orille it 3 dimensions
Table 3 Percentage differences between the experimental freezing and thawing rimes for minced lean beef mulh-dm~cnslonal shapes and the freezing and thawing times calculated by numencal methods Tableau 3 Difference en pourcent entre les temps de congelat ton et de dbcongblatlon experimentaux pour des ]orme~ muh tdtmens umnelles de boeul mulgre hachb et les temps de con4l~lation et de dbcongblation calculbs par des mbthode,s num~rlque~ Data 31
Rectangular brick (4 runs) Two-dimensional irregular shapes (8 runs)
38
Method
Mean 1%)
SD 1%)
Mimmum (",,)
Maximum (",,)
FDM
4.8
90
- 3 ()
17.3
FEMf
12 3
2.9
~2
10.2
Int. J. Refrig. 1987 Vol 10 January
Prediction of freezing and thawing times. D. J. Cleland et al. suggested guidelines were followed. Use of coarse grids, large time steps and/or the simplified FEM formulation did lead to significant inaccuracy in the predictions and therefore worsened overall agreement with experimental freezing and thawing data. The spreads of the FEM predictions were low for all data tested, but mean prediction errors for some threedimensional irregular shapes were offset. This indicated that systematic data errors, probably in the estimation of surface heat transfer coefficients, were important compared w,th random experimental uncertainty.
13
14 15 16
17
Acknowledgements The authors acknowledge the financial support of the Meat Industry Research Institute of New Zealand for this research. In addition, D. J, Cleland acknowledges the support of the New Zealand Meat Producers Board in the form of a Sir Walter Mulholland Fellowship.
18
19
20 21
References 1 2 3
4
5 6
7 8
9
10
11 12
Cleland, A. C., Earle, R. L. Assessment of freezing time predlchon methods J Food S¢t (1984) 49 1034-1042 Cleland, D. J., Cleland, A. C., Earle, R. L. Prediction of freezing and thawing times for foods a review lnt J Re~ly (1986) 9 182 Cleland, A. C.. Earle, R. L. A comparison of analytical and numerical methods for predicting the freezing times of foods J Food Set (1977) 42 1390 1394 Cleland, D. J. Prediction of freezing and thawing times for foods PhD Thems Massey University, New Zealand (1985) Cleland, D. J., Cleland, A. C., Earle, R. L., Byrne, S.J. Prediction of thawing times for food of simple shape lnt J Rel~'iy (1986) 9 220~228 Heldman, D. R. Factors influencing food freezing rates Food Technol (19831 37(4) 103-109 Bonacina, C., Comini, G., Fasano, A., Primieern, M. Numerical solution of phase change problems Int J Heat Mass Transfer (1973) 16 1825-1832 Comini, G., del Guidice, S., Strada, M., Rebeilato, L. The fimte element method in refrigeration engineering lnt J Refrio (1978) I 113-118 Cleland, D. J., Cleland, A. C., Earle, R. L., Byrne, S. J. Predlchon of rates of freezing, thawing and cooling in solids of arbitrary shape using the finite element method lnt J Refri# (1984) 7 6 13 Bonacina, C., Comini, G. On a numerical method for the solution of the unsteady state heat conduction equation with temperature dependent parameters Proc 13th lnt Congr Refriy (1971) 2 329 336 Lees, M. A linear three level difference scheme for quasilinear parabolic equations Math Comput (1966) 20 516--522 Comini, G., del Guidice, S., Lewis, R. W., Zienkiewiez, O. C.
22
23 24
25
26 27
28
29
30
31
Finite element solution of non-hnear heat conduction problems with special reference to phase change lnt J Numer Methods En 9 (1974) 8 613~24 Cleland, A. C., Earle, R. L., Cleland, D. J. The effect of freezing rate on the accuracy of numerical freezing calculations lnt J Relrt9 (1982) 5 294-301 Cleland, A. C., Earle, R. L. The third kind of boundary condition in numerical freezing calculations lnt J Heat Mass Transfer (1977) 20 1029 1034 Comini, G., del Guidiee, S. Thermal aspects of cryosurgery J Heat Trans[br (1976) 98 543-549 Lemmon, E. C. Phase-change techniques of finite element conduction codes in Numeru'al methods m thermal problems (Eds R W. Lewis and K Morgan) Pineridge Press (1979) 149-158 Morgan, K., Lewis, R. W., Zienkiewicz, O. C. An improved algomthm for heat conduction problems with phase change Int J Numer Method.s End1 (1978) 12 (1191 1195 Pham, Q. T. A fast, uncondmonally stable finite-difference scheme for heat conduction with phase change Int J Heat Ma.s,~ Tran.~ler (1985) 28 2079-2084 Pham, Q. T. The use of lumped capacJtance m the finite element solution of heat conduction problems with phase change lnt J Heat Mass Tran,sJer (1986) 29 285 291 Heldman, D. R. Food properties during freezing Food Technol (19821 36(2) 92 96 Cieland, A. C., Earle, R. L. A comparison of methods for predicting the freezing times of cylindrical and spherical foodstuffs J Food Sci (1979) 44 958 963 Cleland, A. C. Heat transfer during freezing of foods and prediction of freezing times PhD Thesl,s Massey UmversJty. New Zealand (1977) Cleland, A. C., Earle, R. L. Pre&ct~on of freezing t~mes for foods in rectangular packages J Food Set (1979) 44 964-970
Frivik, P. E., Thorbergson, E., del Gnidiee, S., Comini, G. Thermal design of pavement structures m seasonal frost areas J Heat TransJer (1977) 99 533 540 Miki, H., Kikukawa, H., Nishimnto, J. Fundamental studies on the thawing of frozen fish 3. Numerical analysis for thawing processes Mere Fac Fi,sherws Kagoshima Umr (1978) 27(1 ) 107115 Hsu, T. R., Pizey, G. On the prediction of fusion rate of ice by finite element analysis J Heat Tran~Jer (1981) 103 727 737
Nishimura, T., Ito, S., Tsuboi, H., Kawamura, Y. Analysis of twodimensional freezing by the finite element method lnt Chem Eng (1985) 25(1) 105 112 Purwadaria, H. K., Heldman, D. R. A finite element model for prediction of freezing rates m food products with anomalous shapes Trans ASAE (1982) 25 827-832 Hayakawa, K., Nonino, C., Sueear, J. Two-dimensional heat conduction m food undergoing freezing: predicting freezing times of rectangular or finitely cylindrical food J Food Scr (1983) 48 1841 1848 Miki, H., Kikukawa, H., Nishimoto, J. Numerical calculations of three-dimensional heat conduction on freezing process of marine foodstuff Bull Jap Soc Set Fisheries (1982) 48(61 775 779
Cleland, D. J., Cleland, A. C., Earle, R. L., Byrne, S. J. Experimental data for freezing and thawing of multi-dimensional objects Int J Refrig (1987) 10 22
Rev. Int. Froid 1987 Vol 10 Janvier
39
An assessment of the accuracy of numerical methods used in the prediction of freezing and thawmg times wa~ made using a comprehensive set of freezing and thawing data for both regular and irregular multidimensional shapes. For regular shapes, a finite difference method gave accurate predictions with reasonable computation costs. Predictions for two finite element method formulations were not always accurate. This was due to practical constraints on the computation costs which meant that time and spatial grids could not always be made sufficiently fine to ensure that the prediction method uncertainty was insignificant compared with the other sources of imprecision. Guidelines are suggested for the use of the finite element method as a freezing or thawing time predictor. These should ensure that the prediction method error is small while keeping the computation costs reasonable.
(Keywords:freezing;thawing;computerprograms)
PrOvision par des m6thodes num riques des temps de congelation et de dbcong61ation pour des formes multidimensionnelles Une Ovaluation de la prbcision de mbthodes numOriques utilisOes pour la provtslon des temps de con qe.lation et de dOeonoblation a btb effeetuOe en utilisant un ensemble complet de donnOes sur la conoblation et la dOcongblatton de formes multidimensionnelles rOguliOres et irrOguliOres. Pour les formes rbguliOres la mOthode des difference~s finies a donnO des provisions prbcises avec des cofits de calcul raisonnables. Les provisions pour deux formules de la mOthode des Olbments finis n'Otaient pas toujours prOcises. C'Otait dfi h des contraintes pratiques sur les cofit~ de calcul, ce qui signifiait que les 9rilles de temps et d'espace ne pouvaient pas toujours ~tre suffisamment fine,s pour s'assurer que rincertitude de la mOthode de provision Otait peu importante en comparaison avec d'autres sources d'impr~cisions. On propose des directives h utiliser dans la mOthode des OlOments finis comme moyen de provision du temps de conoOlation ou de dOcongOlat ion. Avec ces directives l'erreur de la mOthode de provision est faible tout en maintenant des cofits de calcul raisonnables.
(Mots cles: cong~lation; d~cong~latlon; programmes d'ordinateur
Nomenclature Bi D 2D, 3D FDM FEM h k SD
Biot number = hD/k Object thickness or characteristic dimension (m) Two-dimensional, three-dimensional Finite difference method Finite element method Surface heat transfer coefficient ( W m - 2 oC -1 ) Thermal conductivity ( W m -2 '-C 1) Standard deviation of percentage differences
During freezing or thawing of food materials latent heat is released or absorbed over a range of temperatures below O°C. The most physically realistic model for freezing and thawing of foodstuffs is that of transient non-linear heat conduction with temperature-variable thermal properties in the solid object, subject to convective heating or cooling at the object surface1. Without making grossly
* Present address: Auckland University, Auckland, New Zealand 0140-7007/87/010032~)8503.00 Q~ 1987 Butterworth & Co {Publishers} Ltd and IIR
32
Int. d. Refrig. 1987 Vol 10 January
T
Temperature t Ct
Subscripts
a f in x, v, z
Ambient freezing or thawing medium Full unmodified formulation Initial Simplified formulauon In the x, y, or z direction, respectively
simplifying assumptions, the non-linear pamal differential equations representing this model cannot be solved exactly by analytical methods to predict freezing or thawing times. A wide variety of approximate methods have been suggested as solutions 2. These methods fall into two broad groups: numerical methods; and simple formulae. Both groups have advantages and disadvantages so it is expected that food engineers will continue to make use of both for some time to come 1'3 The use of simple formulae is discussed elsewhere 4"5.
Prediction of freezing and thawing times: D. J. Cleland et al. Numerical methods The numerical methods make discrete mathematical approximations to the time and spatial variations defined by the governing partial differential equations for heat conduction. If they are formulated and implemented correctly to reduce numerical truncation and rounding errors, the numerical methods are generally considered the most accurate, reliable and versatile freezing and thawing time prediction methods 1'°. The numerical methods are further subdivided into two major types: finite difference methods; and finite element methods. Both have been extensively used for prediction of freezing or thawing rates and times 2'7'8. Detailed discussion of the advantages and disadvantages of the two methods can be found in previous papers 8"9. The finite difference method (FDM) is limited by practical considerations to use for one-dimensional and regular multi-dimensional shapes where evenly spaced grids can he used, while the finite element method (FEM) can more easily take account of irregular multi-dimensional geometries. Although the spatial discretization is different for the FEM it is usual to employ a FDM scheme for the time progression in transient problems. Bonacina and Comini l° and Cleland and Earle 1 considered different F D M schemes for phase change problems and concluded that the implicit, three time level, Lees scheme 11 was most accurate. Formulations of both the F D M and FEM based on the Lees scheme have been developed and implemented to predict freezing and thawing for a range of geometries s'9' 12.13,14 It is well known that with a very fast computer, unlimited computer memory and appropriate software both time and space intervals can be made sufficiently small that the F D M and FEM will give accurate predictions for multi-dimensional shapes, whether they are irregular or not. But published results indicate that few researchers, let alone industrial food engineers, have such computational resources routinely available. Therefore, it cannot be assumed that the F D M and FEM are universally accurate predictors of freezing and thawing times in practical situations. The size of the available computer memory limits the refinement of the spatial grid that can be used. The sizes of both the time and space intervals are restricted by the computer processing speed and cost. In some situations these limits can be overcome. For example, a F D M used with an 8 × 8 × 8 evenly spaced node grid that required less than 0.5 megabytes of memory and took less than 1000s process time on a Prime 750 computer gave good predictions for freezing of rectangular brick shapes. By comparison, a FEM using a 6 × 6 × 6 evenly spaced nodal grid required about 2 megabytes of memory and 30000s of process time, and yet gave barely adequate prediction accuracy for the same problem 9. Because of such difficulties, modifications to the basic F D M and/or FEM formulations have been developed either to allow the use of a larger time step while maintaining the same prediction accuracy or to reduce the computation effort per time step 9'12't5"1°'17"18'19. With all of these modifications there is a trade-off between the loss or gain in accuracy and the savings or increase in computation costs. In practice, uncertainties in thermal property and boundary conditions data are often substantial, so that provided numerical approximation errors are kept small compared with the physical data
uncertainties, no significant change in prediction accuracy will result. For multi-dimensional problems, users may resort to using coarse spatial grids to save computation time. This can lead to prediction inaccuracy as the spatial approximations may become limiting in accuracy, even when one of the modified methods is used. Another major disadvantage of numerical methods compared with simple formulae is that they require detailed thermal property data over the whole temperature range. There are a number of methods (see, for example, Heldman 2°) to predict thermal properties over a wide range of temperatures from only a few basic measurements, such as moisture content and the initial freezing temperature, so that reasonably accurate prediction of thermal property data is usually possible. Calculation of a heat balance 1'9 as part of the numerical method implementation is one of the best ways to check the accuracy of all the numerical approximations made by the numerical methods 1. The calculated heat flow through the surface of the object should equal the change in the object's internal enthalpy. Close agreement (within ~ 2 °o) should be achieved for a numerical method to be considered both accurate and valid ". Assessment of prediction accuracy o/ numerical methods For problems that can be considered in one dimension (slabs, infinite cylinders and spheres) both the F D M and the FEM are easily implemented, require minimal computational effort and give comparable prediction accuracy. By comparison with experimental data for freezing 3"9'2x and thawing 5 of typical foods, three time level versions of both methods have been proved to be sufficiently accurate predictors of both freezing and thawing times. For freezing of two- and three-dimensional regular shapes (infinite rectangular rods, finite cylinders and rectangular bricks) the accuracy of the Lees scheme F D M has been verified4' 13.22.23. As expected for regular shapes, a full, three time level finite element method formulation (FEMf) was found to give comparable results to the F D M if sufficiently fine time and space grids were used 9. A simplified finite element method formulation (FEMs) 9 was found to be nearly as accurate as the FEMf formulation and gave large savings in computation time. However, it was not clear how the FEM, should be used as, in all cases tested, the spatial grids were reasonably fine and the time steps were controlled so that good heat balances were maintained. For phase change in multi-dimensional irregular shaped objects, comparisons of FEM predictions with experimental data have been made for freezing of roads24, thawing of fish 25, freezing of ice 26'27 and freezing of beefz8'29, all in two dimensions, and for freezing offish 3° in three dimensions. These comparisons were all limited in scope, as they were made for only one or two experimental runs each. Research aims In summary, the accuracy of the three time level F D M has been proved for regular shapes within reasonable computational cost limitations. Although it was expected that the FEM would give accurate prediction of phase change this testing has not been comprehensive,
Rev. Int. Froid 1987 Vol 10 Janvier
33
Prediction of freezing and thawing times: D. J. Cleland et al. Table 1 Percentage differences between the experimental freezing and thawing Umes for Tylose mulU-dlmensumal ,hapes and the freezing and thawing times calculated by numerical methods Tableau 1 Difference en pourcent entre les temps de congelation et de decon.qi'latlon expbrimentaux pour des h~rme~ muhuhmensumnelle~ de l-vlow eL te~ temps de congblation et de dbcongblation calculbs par des mbthodes numbrtque~ Mean (%)
SD (",,1
MmJmum {',,)
Maxtmum ~" ~
Shape
Data 23, 31
Method
Rectangular bricks
Freezing (72 runs)
FDM FEM,
4,4 L9
6.6 7 1
2u 24
Thawing (68 runs)
FDM FEM s
0.6 b3
3 1 34
(, ,~ 12
'~
Freezing and thawing (83 runs)
FEM~ FEM,
0.3 42
5 I 6.1
I! ~ ')~
~l : 2) ~'
Freezmg (42 runs) Thawing (41 runs)
FEMf FEM,
10 26
45 4.8
II ~ ,1
'~, I ? ~)
FEMr FEMs
15 58
55 70
(,
i~ ?1 ,)
Two-&mensmnal irregular shapes
especially for multi-dimensional irregular shapes, and where there are computation cost restrictions. Therefore the aims of this paper are: 1. to assess a three time level finite element method (FEM09 for the prediction of both freezing and thawing times for multi-dimensional irregular shapes by testing against experimental data subject to practical computational cost constraints; 2. to assess the value of the simplified finite element method formulation (FEM,) as proposed by Cleland et al. 9 as an alternative to save computation costs; 3. to formulate guidelines for the use of the FEM for phase change problems that would help to ensure accurate predictions where there are computation cost limitations; and 4. to use data for thawing of rectangular brick shapes to verify further the accuracy of the numerical methods for regular multi-dimensional shapes. Results and discussion
A large set of experimental data was collected for this assessment of the numerical methods. Details of experimental techniques and the experimental results are reported by Cleland et al.31. Assessment of the prediction accuracy of the numerical methods was made by direct comparison of the numerically calculated times for the thermodynamic centre temperature of the objects to reach - 1 0 ° C for freezing, or 0°C for thawing, with the experimentally measured times. The percentage differences were calculated from: 3/0 Difference =
Predicted time-Experimental time 100 ×--Experimental time 1 (1)
In all cases the thermal property data used for the numerical method predictions were those given by Cleland et al. s and at least 1000 time steps were used in each numerical calculation. Prediction of freezing and thawing times for rectangular bricks Table I summarizes the prediction accuracy of the
34
Int. J. Refrig. 1987 Vol 10 January
,~ s
-, ,
F D M 23 and the simphfied finite element method formulation, F E M f , compared with experimental data for freezing23 and thawing 31 of rectangular brick shapes of Tylose. The full finite element method formulation. FEMf 9, was not used for the full data set because computation costs were large with the grids used. If the FEMf had been used, the results would be expected to match the F D M results 9 closely. Because of symmetry only an octant of each rectangular brick was modelled. Within this octant an evenly spaced 9 x 9 x 9 node grid was used in the F D M calculations. The FEM calculations used a 7 ~ 7 x 7 node grid with 216 evenly sized, eight-node, linear isoparametric, brickshaped elements. A finer grid was not practicable due to limitations on computation cost and computer memory size. For all except a few runs the calculated heat balance/'9 agreed to within 2 o . Computer process times on a Prime 750 computer were typically of the order of 1000s per run for the F D M and 10000s per run for the FEM,. The FEMr was about five times slower than the FEM,. The freezing results have been discussed previously ". The high spread of percentage differences and the overall under-prediction was attributed to freezing rate affecting the precision of the thermal property data rather than to uncertainty in the experimental data or imprecision due to the numerical methods 13. Thawing is inherently slower than freezing, so the same effect was not expected for the thawing data. For thawing, Table 1 shows that the FEMs formulation had only a slightly greater offset of the mean prediction error than the predictions for the FDM. As was the case for thawing of the simple slab, infinite cylinder and sphere shapes 5, this difference was thought to be due to the relatively crude way in which variations in the thermal properties are incorporated into the FEMs. The relatively coarse nodal grid used for the three-dimensional geometry accentuated this problem, so the difference between the F D M and FEM, is greater for rectangular bricks than for the simple one-dimensional shapes. Predicted temperature profiles agreed closely with experimentally measured profiles, as illustrated in Figure 1. Uncertainty in thermocouple placement was probably a major contributor to the differences between them. Both the F D M and FEM, formulations predicted oscillating
Prediction of freezing and thawing times. D. J. Cleland et al. 40 I
I E
0
-20 0
I
I
I
I
1
2
5
4
Time (hrs) Figure I Typical temperature l,ersu.~ time profiles for thawing a rectangular brick of Tylose. Run 7 (Table 1 of Reference 31) D~= 0.104m; D,.=0.125m: D~=0.152m; h = 4 1 . 0 W m -2 'C - i , Ta= 45 0 C; T~. = - 12.0 C. , Experimentally measured: . predicted by F D M : . predicted by F E M s. 1. Temperature at the centre of the :,,z surface: 2. temperature at the geometric centre Figure 1 Dtagramme type tempbrature temps pour la dbcong~latton d'un bloc parallel epipbdtque de 7].lose. Essat 7 (le tableau I de la r~ference 31). Dx=O,lO4m. D y - O , 1 2 5 m . Dz=O,152m; h = 4 1 . O W m 2C 1. T , = 4 5 . 0 C. T,, = - 12.0 C --. Me~ur~ e,:p~rimentalement. . prbt u par F D M. . prbt,u par F E M , 1. Tempbrature au centre de la ~urface x.z. 2, temperature au centre geomi, trique
surface temperatures due to the discrete nature of the numerical methods. For the FEM~ formulation this was magnified by the coarse grid used for the predictions 4. Use of the three time level enthalpy method ~8 or the other modifications to the method to improve estimation of the thermal properties in the phase change temperature range ~2'15 ~ may help reduce this type of oscillation. The differences between predicted and measured thawing times were not statistically correlated with any independent variable in the experiments. Further, the 95 °'o confidence bounds for the thawing time predictions of the numerical methods are similar in magnitude to the estimated experimental uncertainty bounds of + 5.4 °,o31. This suggests that for the F D M and FEMr formulation, most of the difference between the predictions and experimental values was due to experimental uncertainty and not imprecision resulting from the implementation of the numerical methods themselves. The only significant inaccuracy arising from the numerical methods was that introduced by the crude thermal property estimation in the FEM, formulation, thus indicating that care should be exercised in using the simplified formulation with coarse grids.
Predictions fi~r two-dimensional irregular shapes The experimental freezing and thawing times for twodimensional irregular Tylose shapes are given in Table 2 of Reference 31. The grids shown in Figures 1-7 of Reference 31 were used to make predictions by both FEM formulations. Eight-node, quadratic, isoparametric, curved rectangular elements were used. The grids prepared ensured that at least five elements (11 nodes) existed between the thermodynamic centre and the outside surface of each object. The FEMf formulation
gave satisfactory heat balances for all the runs, but the FEM s formulation did not always achieve heat balances within the 2 ° o error bounds, although computation costs were reduced by about 80°4, for the same time step size. A typical calculation for the FEMf required about 2 megabytes of memory and took about 50000s process time on a Prime 750 computer. The results, summarized in Table 1, show that both FEM formulations gave accurate predictions of freezing and thawing times (within the tolerances allowed for experimental error) for all the shapes tested. The 95 ° o confidence bounds for the percentage differences of the FEMf formulation of - 9 . 7 to 10.3 ° o are similar to the estimated experimental error of + 8.9 ",,31 indicating that most of the difference was explainable by the experimental uncertainty. Again, the FEM, formulation did introduce some extra inaccuracy due to the cruder thermal property estimation, as indicated by the slightly higher spread and mean offset of the percentage differences. The freezing data were under-predicted compared with the thawing data by both formulations. This was attributed to the freezing rate effect t3 rather than numerical approximation problems because good heat balances were consistently achieved for the FEMf formulation. Typical predicted temperature profiles (Figures 2 and 3) compared well with those measured experimentally. Again, differences between them are probably due in part to uncertainty in thermocouple placement. Also, the possibility of undesired variations m the local surface heat transfer coefficient for these experiments 4'5 meant that measured temperatures at points on, or near, the surface had significant uncertainties associated with them. The surface temperatures predicted by the FEMs formulation
,0 7 20-
oQ_
E -20 i
- 4 0 ~ - -
0
- - F
2
-
-
T . . . . . .
I
4
6
--
---q
8
Time (hrs) Figure 2
Typical temperature rersu~ time profiles for freezing a Tylose two-dimensional irregular shape. Run 29 (Table 2 of Reference 31), s h a p e n o 3. D = 0 . 1 3 7 m ; h = 2 8 . 0 W m - 2 C 1; T a = _ 3 0 . 5 o C , Tm= 32.8 C. - - , Experimentally measured: , predicted by FEMf, -- , predicted by F E M s. The object shape and thermocouple positions 2, 3 and 5 are shown in Figure 3 of Reference 31 Figure 2 Dlagramme type tempbrature temp,s pour la congblation d'une forme irr#guh#re de Tylose gt deux dimensions. Essal 29 (tableau 2 de rbference 31). Forme no. 3 D = O , 1 3 7 m . h = 2 8 , 0 W m 2 'C I T, = - 30.5'~C. T~ = 32,8 'C - - , Mesur~ expbrimentalement. -. prbvu par F E M f: -. prei,u par F E M s. La forme de l'objet et les emplacements du thermocouple sont indtqukes .figure 3 de la rbference 31
Rev. Int. Froid 1987 Vol 10 Janvier
35
Prediction of freezing and thawing times." D. J. Cleland et al. heat balances were not achieved for all of the predictions. mainly because spatial grids were coarse for all of the calculations with three-dimensional grids. The size of the available computer memory, the cost of the computer processing and/or the large data preparation time meant than finer grids could not be used. The grid used for the pyramid shape is shown m Figure 8 of Reference 31. It is basically the same 7 x 7 × 7 node grid used for the rectangular brick predictions, but is distorted to the pyramid dimensions. Because of the three-dimensional irregularity a full quadrant of the shape was modelled. Consequently the grid was coarse. with only three linear elements between surface and centre along the height of the pyramid. Computation costs were the same as for the rectangular brick predictions. The predictions have similar spread to those for freezing and thawing of rectangular bricks but the mean was offset by about 10°o . The low standard deviation of the prediction errors suggests that the source of this error is systematic, not random. Two possible reasons are: the coarse grid leads to significant prediction inaccuracy: and the assumption that h was the same for the pyramid and the reactangular brick shapes 3~ may not hold. Both reasons probably contributed. In the numerical calculations, freezing and thawing were considered complete when the slowest cooling or heating node reached the desired final temperature. If the true thermodynamic centre of the object does not exactly coincide with the slowest node position at this time, then the object will not be completely frozen or thawed. Under-prediction of the time for the complete process will result from this effect. The effect is smallest for fine finite element grids, but for coarse, poorly positioned grids the extra time could be substantial. Further, linear elements magnify the problem because they use linear approximation functions between nodes. The true temperature profile and any local maximum or minimum values of the temperature within each element may not be properly modelled. Quadratic and higher-order elements allow more accurate approximations to the temperature profile to be made and therefore often give more accurate
oscillated significantly compared with the predictions for the FEMf formulation. Again, this difference was attributed to the crude thermal property estimation used in the simple formulation, which was accentuated by the relatively small number of quadratic elements in the grids used for the predictions.
Predictions for three-dimensional irregular shapes Table 2 summarizes the percentage difference between experimental freezing and thawing times for some threedimensional irregular shapes (Table 3 of Reference 31) and the predictions by the two FEM formulations. Good
10
? v t-
Ill
t:l.
E
-10-
-20 --~ I
-30-~--0
I
[. . . . .
2
4
I
I
I
6
8
10
.....
7
12
Time (hrs) Figure 3 Typical t e m p e r a t u r e versus t i m e profiles for t h a w i n g a Tylose t w o - d i m e n s i o n a l irregular shape. Run 31 (Table 2 of Reference 31 ), s h a p e no. 4: D = 0 . 0 9 9 m : h = 2 8 . 0 W m - 2 " C - ~ : T a = 2 0 . 1 ° C : Tin= - 21.5'-C. - - , Experimentally measured: , predicted by F E M f , -, predicted by F E M s. The object shape a n d t h e r m o c o u p l e positions 2, 3 a n d 6 are s h o w n in F i g u r e 4 of Reference 31 F i g u r e 3 Dlagramme type tempbrature temps pour la dbconoblatum
d'une forme irrboulibre de Tylose h 2 dimensions Essat 31 (tableau 2 de rbference 31), forme no. 4 D=O,O99m. h = 2 8 , 0 W m -2 C -1 T = 2 0 , 1 ~ C . T~ = - - ~1,,~ C , Mesurb experirnentulement, prbvu par F E M y. -- -, prbvu par FEM,. La forme de l'objet et les emplacements du thermocouple sont indlqubs figure 4 de rbfferenee 31
Table 2 Percentage differences between the e x p e r i m e n t a l freezing and t h a w i n g u m e s for l y l o s e t h r e e - d i m e n s i o n a l irregular shapes and the freezing a n d t h a w i n g times c a l c u l a t e d by n u m e r i c a l m e t h o d s T a b l e a u 2 Difference en pourcent entre les temp~ de congblatlon et de dbcon,qielatton expermwntaux pour des ]brines muh Mlmen~ionnelles lrrequhere~s dc
Tylose et les temps de congblation et de dbcongblation calculbs par des mbthodes numbriques
D a t a 31
Grid
Method
Mean I ",,)
Pyramid {6 runsl
3D
FEMf FEM,
-- 10.3 - 14 5
6.7 4 8
IS 2 _t1,_ ~ ~
- 20 8q
3D
FEMf FEM~
- 5.5 --I 2
6 8 12 7
- t2 I lt~ I
5.3 13 9
2D
FEMf FEM s
0.6 7.7
70 8.3
,~ 7 J6
13.9 18 2
ID
FEMf FEM s
0.4 2,6
67 65
,' v 4.9
9 6 11 6
3D
FEMf FEM s
5.2 9 8
13 0 14 9
-- 12.9 - ~) 4
25.1 28.8
2D
FEMf FEMs
3,7 16.4
13.8 18 7
- lh 1 4 1
24 7 48 4
3D
FEMf FEM s
35.3 33.5
89 16.3
20 8 8.6
44 2 58.4
Sphere (6 runs)
Egg
(6 runs)
Fish (6 runs)
36
Int. J. Refrig. 1987 Vol 10 January
SD I",)
Minimum 1" ~
Maximum I ",,t
Prediction of freezing and thawing times. D. J. Cleland et al. predictions. FEM calculations for the pyramid shape were affected to some extent in the manner just described as a reasonably coarse grid of linear elements was used, The estimate of the position of the thermodynamic centre was not accurate so the 'centre' node position was not very close to the actual thermodynamic centre. The amount of offset was greater for the freezing runs where the Biot number, Bi, was low, suggesting, according to the criteria of Cleland and Earle ~, that the second reason (systematic error in h) was also an important contributor. The other three shapes were a sphere, an egg and a fish shape. The sphere is a regular shape and can be modelled by one- or two-dimensional grids rather than threedimensionally if axisymmetrical formulations of the FE M are used. Three different grids were used (Figures 1 and 11 of Reference 31 i: 1. an 11 node, 10 linear element, one-dimensional grid (typical computer process time of 100 s per run on a Prime 750 computer): 2. a two-dimensional grid modelling a quadrant of the cross section (a circle) with 21, eight-node, quadratic, isoparametric elements (typical process time of 5000 s): and 3. a 208 node, three-dimensional grid modelling an octant of the sphere with 27 quadratic, isoparametric elements (typical process time of 70 000 s), The three-dimensional grid was coarse, with only three quadratic elements between the surface and the centre. These were arranged in three expanding shells so that the elements were all approximately the same size and none was too distorted from the basic rectangular brick shape. The sphere shape predictions best illustrate the effect of using a coarse grid and the difference between the FEMf and FEM~ formulations. The results generated using both the finer one- and two-dimensional grids with the FEMf formulation were very similar. The results generated using the coarser three-dimensional grid were significantly less accurate. Incorrect implementation of the three-dimensional part of the program was unlikely to be the cause as the program had been extensively tested 9. The position of the thermodynamic centre was known exactly for the sphere, and quadratic elements were used. Therefore, the problem with the coarse grid was less acute than for the pyramid shape predictions and mean underprediction for the coarser three-dimensional grid, compared with the finer one- and two-dimensional grids, was less than the mean under-prediction for the pyramid shape. For the sphere, the FEMs formulation gave comparable results with the formulation, but where the element size was larger in the two- and threedimensional grids the results had a larger variation due to the less representative incorporation of thermal properties into the finite element scheme. For both finite element methods the freezing and thawing time prediction accuracy as well as the prediction of temperature profiles was similar to that achieved for a large set of sphere thawing data 5. The experimental techniques were similar, except for use of the plastic moulding for the shapes, so similar prediction accuracy was expected. The egg shape has one axis of rotational symmetry so it could be modelled by either a two- or three-dimensional grid (Figure 9 of Reference 31). Both the grids were
double versions of the sphere grids, with one half elongated to fit the ovoid shape. The fish shape, although irregular in all three dimensions, was closely approximated by a shape with two planes of symmetry, so only a quadrant was modelled. The grid is shown in Figures l0 and 12 of Reference 31. It is the same basic three-dimensional grid used for the egg shape, but the nodal coordinates, especially for the surface nodes, were modified to fit the actual shape. A typical calculation using the threedimensional egg grid or for the fish shape took 150000s process time on a Prime 750 computer and required 4 megabytes of memory. For the egg and fish shapes it was assumed that the surface heat transfer coefficient was the same as that measured on the moulded plastic spheres 31. The additional inaccuracy introduced by this assumption could not be easily assessed. This extra experimental uncertainty, plus the uncertainty due to difficulties in measuring and describing these shapes mathematically, as well as the approximation of the shape with a coarse FEM grid, meant that poorer agreement was expected for the egg and fish shapes than for the sphere and pyramid shapes. For the egg shape the mean prediction error was close to zero, so the assumption that the surface heat transfer coefficient was the same as it was for the sphere seemed reasonable. The higher standard deviation of the predictions for the egg shape compared with the pyramid and sphere shapes was considered to be mainly the result of the additional uncertainty in the experiments rather than uncertainty introduced in the implementation of the FEM. Ifa finer grid had been used, the FEMf would not be expected to predict the experimental data for the egg much better. The accuracy of temperature profiles for the egg shape (Figure 4), calculated using the three-dimensional grid reflected the accuracy of the freezing and thawing time predictions. Temperature profiles for some surface positions may have been affected by localized variation in h 4. The predicted temperature profiles were affected by the coarse three-dimensional grids used, and therefore displayed some oscillatory behaviour. As expected, this was particularly obvious for the FEM, formulation predictions, for positions near the object surface and at temperatures close to the phase change temperature range 9. For the fish shape the predictions of freezing and thawing times were consistently high by 20-40 °,o but the standard deviation of prediction errors was comparable with that for the egg shape. The offset was probably caused by systematic error in the surface heat transfer coefficient, which may not have been the same as that for the sphere or egg due to: thinner plastic thickness arising during moulding; and/or better thermal contact between the Tylose and the plastic mould for the fish shape than for the sphere and egg shapes. As discussed previously 31, these factors could not be quantitatively assessed as it was not possible to measure h directly for the fish shape. However, measurement of the plastic thickness showed variations of up to 15 O/o of the average total thickness over the surface of the fish shape and differences of up to 10 °/o in the thickness between the fish and sphere shapes. Further, predicted temperature profiles for the fish shape 4 showed trends consistent with a systematic underprediction of h.
Rev. Int. Froid 1987 Vol 10 Janvier
37
Prediction of freezing and thawing times. D. J. Cleland et al. and the two-dimensional irregular shapes. This was attributed to systematic error in the estimate of the surface heat transfer coefficient. The assumption that h was the same for the meat experiments as for the Tylose experiments may not have been valid 3~. Tylose and minced lean beef have different textures, so the contact resistance between the object mould walls and the materials may not have been equal. In view of the small data set, firm conclusions cannot be made about prediction accuracy for the meat data, Howe~er, the general agreement should be noted
Prediction of freezing and thawing data for minced lean beef Freezing and thawing data for minced lean beef rectangular bricks and two-dimensional irregular shapes (Table 4 of Reference 31) were also predicted using the FDM and FEMf. The same grids employed previously were used for the calculations. Table 3 summarizes the results. It was expected that the numerical methods would introduce negligible additional prediction inaccuracy when applied to freezing and thawing of actual foods. However, an increase in thermal property data and experimental uncertainty relative to the Tylose data was expected4. The larger standard deviation of the prediction errors for the bricks compared with the equivalent Tylose results may reflect the difficulties in experimenting with actual food materials. The mean percentage differences for the meat data were higher than for Tylose by about 10% for both the bricks
Guidelines ,[br use of the finite element method Some guidelines for application of the FEM were established by Cleland et al. 9. In the light of the present assessment against experimental data these can be refined. Adherence to these guidelines should lead to predictions in which the error arising from the numerical approximations is negligible compared with experimental and data uncertainty while keeping computation costs as low as possible. The guidelines are'
30-
1. to use sufficiently small nine steps that the heat balance agrees to within 2".~: 2. to use at least 8-10 nodes and at least five elements between the thermodynamic centre and the surface in the FEM grids; 3. to position enough nodes near the estimated position of the thermodynamic centre to ensure that one will ultimately be reasonably close; and 4. to use the full FEM formulation if computational cost limitations allow.
20-.
o
~
-10 -
Of these guidelines, the first two are the most important. In almost all situations, using more refined grids and smaller time steps than those suggested will only give small increases in accuracy and will probably not justify the increase in computation costs
-20 -
-300
1
I
1
I
1
1
2
3
4
5
I. . . .
q
6
7
Conclusions
Time (hrs)
Practical constraints on computation power mean that users of numerical methods, especially the FEM, will need to choose the grid size carefully. The finite difference method used for regular two- and three-dimensional shapes predicted thawing times a~ accurately as it predicted freezing times. Even with computation cost constraints the method could be implemented with negligible inaccuracy arising from the numerical approximations used. The implementations of the FEM used for both regular and irregular multi-dimensional shapes did not lead to significant inaccuracy in predictions provided the
Figure 4 Typical temperature ver,~us time profiles for freezing or thawing a Tylose three-dimensional irregular shape. Run 16 [Table 3 of Reference 31), egg shape: D = 0 A 7 0 m ; h = 5 1 . 4 W m - 2 C ~, T a = - 2 2 . 6 ° C ; Tm=20.6°C. - - , Experimentally measured. . predicted by FEMf; - - - - - , predicted by F E M s. The object shape and thermocouple positions 2, 4 and 6 are shown in Figure 4 of Reference 31. calculated using the three-dimensional grid Figure 4 Diagramme type tempbrature-temps pour la congblatton et la dbcong~lation d'une forme irr~gulibre de Tylose it 3 dimensions. Essai 16 (tableau 3 de rbference 31). Forme de l'oeuf D=O,17Om. h = 5 1 , 4 W m - 2 o c 7, T = _ 2 2 , 6 o c . T =20,6~C. , Me,surb expbrimentalement, - - , prbvu par F E M y, , prbvu par F E M s La forme de l' objet et les emplacements du thermoeouple sont indiqubs figure 4 de la rbference 31, calculbes en utilisant la orille it 3 dimensions
Table 3 Percentage differences between the experimental freezing and thawing rimes for minced lean beef mulh-dm~cnslonal shapes and the freezing and thawing times calculated by numencal methods Tableau 3 Difference en pourcent entre les temps de congelat ton et de dbcongblatlon experimentaux pour des ]orme~ muh tdtmens umnelles de boeul mulgre hachb et les temps de con4l~lation et de dbcongblation calculbs par des mbthode,s num~rlque~ Data 31
Rectangular brick (4 runs) Two-dimensional irregular shapes (8 runs)
38
Method
Mean 1%)
SD 1%)
Mimmum (",,)
Maximum (",,)
FDM
4.8
90
- 3 ()
17.3
FEMf
12 3
2.9
~2
10.2
Int. J. Refrig. 1987 Vol 10 January
Prediction of freezing and thawing times. D. J. Cleland et al. suggested guidelines were followed. Use of coarse grids, large time steps and/or the simplified FEM formulation did lead to significant inaccuracy in the predictions and therefore worsened overall agreement with experimental freezing and thawing data. The spreads of the FEM predictions were low for all data tested, but mean prediction errors for some threedimensional irregular shapes were offset. This indicated that systematic data errors, probably in the estimation of surface heat transfer coefficients, were important compared w,th random experimental uncertainty.
13
14 15 16
17
Acknowledgements The authors acknowledge the financial support of the Meat Industry Research Institute of New Zealand for this research. In addition, D. J, Cleland acknowledges the support of the New Zealand Meat Producers Board in the form of a Sir Walter Mulholland Fellowship.
18
19
20 21
References 1 2 3
4
5 6
7 8
9
10
11 12
Cleland, A. C., Earle, R. L. Assessment of freezing time predlchon methods J Food S¢t (1984) 49 1034-1042 Cleland, D. J., Cleland, A. C., Earle, R. L. Prediction of freezing and thawing times for foods a review lnt J Re~ly (1986) 9 182 Cleland, A. C.. Earle, R. L. A comparison of analytical and numerical methods for predicting the freezing times of foods J Food Set (1977) 42 1390 1394 Cleland, D. J. Prediction of freezing and thawing times for foods PhD Thems Massey University, New Zealand (1985) Cleland, D. J., Cleland, A. C., Earle, R. L., Byrne, S.J. Prediction of thawing times for food of simple shape lnt J Rel~'iy (1986) 9 220~228 Heldman, D. R. Factors influencing food freezing rates Food Technol (19831 37(4) 103-109 Bonacina, C., Comini, G., Fasano, A., Primieern, M. Numerical solution of phase change problems Int J Heat Mass Transfer (1973) 16 1825-1832 Comini, G., del Guidice, S., Strada, M., Rebeilato, L. The fimte element method in refrigeration engineering lnt J Refrio (1978) I 113-118 Cleland, D. J., Cleland, A. C., Earle, R. L., Byrne, S. J. Predlchon of rates of freezing, thawing and cooling in solids of arbitrary shape using the finite element method lnt J Refri# (1984) 7 6 13 Bonacina, C., Comini, G. On a numerical method for the solution of the unsteady state heat conduction equation with temperature dependent parameters Proc 13th lnt Congr Refriy (1971) 2 329 336 Lees, M. A linear three level difference scheme for quasilinear parabolic equations Math Comput (1966) 20 516--522 Comini, G., del Guidice, S., Lewis, R. W., Zienkiewiez, O. C.
22
23 24
25
26 27
28
29
30
31
Finite element solution of non-hnear heat conduction problems with special reference to phase change lnt J Numer Methods En 9 (1974) 8 613~24 Cleland, A. C., Earle, R. L., Cleland, D. J. The effect of freezing rate on the accuracy of numerical freezing calculations lnt J Relrt9 (1982) 5 294-301 Cleland, A. C., Earle, R. L. The third kind of boundary condition in numerical freezing calculations lnt J Heat Mass Transfer (1977) 20 1029 1034 Comini, G., del Guidiee, S. Thermal aspects of cryosurgery J Heat Trans[br (1976) 98 543-549 Lemmon, E. C. Phase-change techniques of finite element conduction codes in Numeru'al methods m thermal problems (Eds R W. Lewis and K Morgan) Pineridge Press (1979) 149-158 Morgan, K., Lewis, R. W., Zienkiewicz, O. C. An improved algomthm for heat conduction problems with phase change Int J Numer Method.s End1 (1978) 12 (1191 1195 Pham, Q. T. A fast, uncondmonally stable finite-difference scheme for heat conduction with phase change Int J Heat Ma.s,~ Tran.~ler (1985) 28 2079-2084 Pham, Q. T. The use of lumped capacJtance m the finite element solution of heat conduction problems with phase change lnt J Heat Mass Tran,sJer (1986) 29 285 291 Heldman, D. R. Food properties during freezing Food Technol (19821 36(2) 92 96 Cieland, A. C., Earle, R. L. A comparison of methods for predicting the freezing times of cylindrical and spherical foodstuffs J Food Sci (1979) 44 958 963 Cleland, A. C. Heat transfer during freezing of foods and prediction of freezing times PhD Thesl,s Massey UmversJty. New Zealand (1977) Cleland, A. C., Earle, R. L. Pre&ct~on of freezing t~mes for foods in rectangular packages J Food Set (1979) 44 964-970
Frivik, P. E., Thorbergson, E., del Gnidiee, S., Comini, G. Thermal design of pavement structures m seasonal frost areas J Heat TransJer (1977) 99 533 540 Miki, H., Kikukawa, H., Nishimnto, J. Fundamental studies on the thawing of frozen fish 3. Numerical analysis for thawing processes Mere Fac Fi,sherws Kagoshima Umr (1978) 27(1 ) 107115 Hsu, T. R., Pizey, G. On the prediction of fusion rate of ice by finite element analysis J Heat Tran~Jer (1981) 103 727 737
Nishimura, T., Ito, S., Tsuboi, H., Kawamura, Y. Analysis of twodimensional freezing by the finite element method lnt Chem Eng (1985) 25(1) 105 112 Purwadaria, H. K., Heldman, D. R. A finite element model for prediction of freezing rates m food products with anomalous shapes Trans ASAE (1982) 25 827-832 Hayakawa, K., Nonino, C., Sueear, J. Two-dimensional heat conduction m food undergoing freezing: predicting freezing times of rectangular or finitely cylindrical food J Food Scr (1983) 48 1841 1848 Miki, H., Kikukawa, H., Nishimoto, J. Numerical calculations of three-dimensional heat conduction on freezing process of marine foodstuff Bull Jap Soc Set Fisheries (1982) 48(61 775 779
Cleland, D. J., Cleland, A. C., Earle, R. L., Byrne, S. J. Experimental data for freezing and thawing of multi-dimensional objects Int J Refrig (1987) 10 22
Rev. Int. Froid 1987 Vol 10 Janvier
39