Tính toán thời gian cấp đông thực phẩm dạng trụ vô hạn và cầu

ISSN 1859-1531 - TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ ĐẠI HỌC ĐÀ NẴNG, SỐ 11(96).2015, QUYỂN 2 11 TÍNH TOÁN THỜI GIAN CẤP ĐÔNG THỰC PHẨM DẠNG TRỤ VÔ HẠN VÀ CẦU CALCULATION FREEZING TIME FOR INFINITE CYLINDER AND SPHERES SHAPED FOOD Nguyễn Bốn1, Võ Chí Chính1, Hoàng Minh Tuấn2 1Trường Đại học Bách khoa, Đại học Đà Nẵng; vcchinh@dut.udn.vn 2Trường Cao đẳng Công nghiệp Huế; hmtuan@hueic.edu.vn Tóm tắt - Hiện nay, điện năng tiêu tốn cho các thiết bị lạnh đông chiếm trên 38% điện năng tiêu thụ trong các nhà máy chế biến thủy sản. Do đó, việc dự báo chính xác thời gian cấp đông thực phẩm có ý nghĩa quan trọng trong việc tiết kiệm năng lượng, giảm chi phí sản xuất, nâng cao chất lượng và tăng tính cạnh tranh của thực phẩm đông lạnh. Một phương pháp đơn giản để tính dự đoán thời gian cấp đông của thực phẩm dạng trụ vô hạn và cầu được đề xuất trong bài báo này. Phương pháp này dựa vào phương trình cân bằng nhiệt tức thời của vật cho các giai đoạn chuyển pha với quá trình truyền nhiệt không ổn định trong giai đoạn làm lạnh, chuyển pha và quá lạnh. Phương pháp này cho kết quả khá chính xác so với các phương pháp giải tích đã có từ trước đến nay. Tất cả các tính toán này dễ dàng được lập trình trên máy tính và dự đoán chính xác thời gian đông lạnh thực phẩm. Abstract - Nowadays, electrical power using for freezing equipment accounts for over 38% of total power of seafood processing factories. Therefore, the exact prediction of food freezing time plays an important role in saving energy, reducing production cost, improving quality and increasing competition of food freezing. A simple method used to calculate freezing time of infinite cylinders and spheres shape food was proposed in this article. This model is based on the energy balance equation of food products for transition phase with unsteady state heat transfer solutions in pre-cooling, phase change and tempering time. This method gave more accurate results than the previous methods, including the finite difference method and finite element. All calculations are easily programmable on computer and predict freezing time of food accurately. Từ khóa - trụ vô hạn; cầu; thời gian đóng băng; thực phẩm; truyền nhiệt không ổn định. Key words - infinite cylinders; spheres; freezing time; food product; transient heat transfer. 1. Đặt vấn đề Nhiều phương pháp giải tích đã được đề xuất để dự đoán thời gian đóng băng cho thực phẩm có hình dạng đơn giản. Hầu hết các phương pháp này được phát triển từ phương trình của Plank (1941). Những hạn chế trong phương trình Plank là nhiệt độ ban đầu của vật chính bằng nhiệt độ điểm đông, tính chất nhiệt vật lý không thay đổi, dẫn nhiệt trong lớp băng là ổn định, không có hai giai đoạn làm lạnh và quá lạnh. Các phương pháp khác của Nagaoka và cộng sự (1955); Plank (1963); Cleland và Earle (1977, 1982); Hung và Thompson (1983) đã được đề xuất để điều chỉnh phương trình Plank (1941). Bên cạnh đó, một số phương pháp của Mascheroni và Calvel (1982); De Michelis và Calveo (1982, 1983); Castaigne (1985a, b); Castaigne và Lacroix (1985) có kết hợp với phương trình cân bằng nhiệt trong quá trình làm lạnh và quá lạnh. Những phương pháp trên thường được dùng để dự đoán cho vật ẩm có hình dạng tấm phẳng. Dự đoán thời gian cấp đông cho vật có dạng hình trụ vô hạn hoặc hình cầu là tỷ số thời gian cấp đông của vật có dạng tấm phẳng trên hệ số quy đổi thời gian đông lạnh của vật hay còn gọi là kích thước truyền nhiệt tương đương. Hệ số này được phát triển bởi Cleland và cộng sự (1987a, b); McNabb (1990b); Hossain (1992b). Việc sử dụng hệ số quy đổi thời gian đông lạnh rất dễ dàng và nhanh chóng để dự đoán thời gian cấp đông cho các vật có hình dạng khác. Tuy nhiên, nó có thể mắc phải sai số trong quá trình quy đổi. Trong khuôn khổ bài báo này, tác giả đề xuất phương pháp dự đoán thời gian đóng băng cho thực phẩm có hình dạng trụ vô hạn và cầu mà không dùng hệ số quy đổi. 2. Kết quả nghiên cứu và khảo sát 2.1. Mô tả quá trình đông lạnh thực phẩm Quá trình đông lạnh thực phẩm được mô tả bởi đường đặc tính t(t) theo thời gian t như Hình 1, gồm ba giai đoạn phân biệt: Hình 1. Đặc tính t(t) của quá trình cấp đông thực phẩm 1. Hạ nhiệt độ của thực phẩm từ nhiệt độ ban đầu ti đến nhiệt độ bắt đầu đóng băng t0; 2. Chuyển pha từ lỏng sang rắn của phần ẩm trong thực phẩm, ở nhiệt độ t0 = const, tỏa nhiệt rc; 3. Quá lạnh thực phẩm ở pha rắn. Lúc này, thực phẩm được tiếp tục làm lạnh từ nhiệt độ t0 đến khi nhiệt độ tâm vật đạt tc theo yêu cầu. Thời gian cấp đông t là tổng thời gian của ba giai đoạn trên 1 2 3t t t t   , [s]. Trong công nghiệp thực phẩm, việc giảm t , nhất là 2t , sẽ làm giảm tiêu hao năng lượng khi cấp đông và thường làm tăng chất lượng của thực phẩm đông lạnh. 12 Nguyễn Bốn, Võ Chí Chính, Hoàng Minh Tuấn 2.2. Phát biểu bài toán đông lạnh thực phẩm và các giả thuyết nghiên cứu 2.2.1. Phát biểu bài toán Xét thực phẩm dạng trụ vô hạn hoặc cầu có kích thước như Hình 2. Khối lượng riêng và nhiệt dung riêng pha ẩm và pha rắn là l , lc và s , sc . Thực phẩm có độ ẩm  , hệ số dẫn nhiệt l và s , nhiệt độ ban đầu it , nhiệt độ bắt đầu đóng băng 0t nhiệt hóa rắn pha ẩm rc và được quá lạnh để nhiệt độ tâm của vật đạt ct bằng cách cho vật tiếp xúc với môi trường có nhiệt độ 0f c it t t t   với hệ số tỏa nhiệt . Cần tính thời gian cấp đông t theo các thông số của bài toán. Hình 2. Trường  ,t r t trong vật ẩm dạng trụ vô hạn 2.2.2. Giả thiết nghiên cứu 1. Tại mỗi thời điểmt coi nhiệt độ t(t) và các thông số vật lý , ,c  của vật ẩm là phân bố đều trong vật. 2. Quá trình đóng băng là quá trình chuyển pha, từ lỏng sang rắn của các thành phần ẩm  trong thực phẩm, ở nhiệt độ bắt đầu đóng băng t0 = const và nhiệt hóa rắn pha ẩm rc = const, lấy trung bình cho các thành phần ẩm của thực phẩm. 2.3. Tính toán thời gian cấp đông vật ẩm dạng trụ vô hạn 2.3.1. Tính thời gian làm lạnh vật ẩm 1t Thời gian làm lạnh vật ẩm 1t từ nhiệt độ ban đầu ti đến nhiệt độ hóa rắn t0 được tính theo phương trình cân bằng nhiệt tức thời lúc t cho 2. .V R l sau thời gian vô cùng bé  d dt t t t   khi  t t d tt   là:    . . . . . .l l fV c dt t t F d  t   (1)   0 . . . . i t l lft dt F d V ct t t  t           . . . .. l l F V c f i ft t t t e  t t      (2)  1 0 . . .ln , . i fl l f t tV c s F t t  t      (3) 2.3.2. Tính thời gian chuyển pha 2t Gọi r là bán kính lớp băng đã tạo ra trước thời điểm t , dr là lớp băng mới tạo ra sau thời gian dt . Phương trình cân bằng nhiệt tức thời lúc t cho khối băng hình trụ mới tạo ra dV = F.dr = 2.π.r.l.dr sau thời gian vô cùng bé dt là:  0 2. . . . . . . . 1 1 .ln 2. . 2. . . f l c l t t r l dr r ld R r R    t         (4)  0 . . . ln ln . . .. l c l l f r d R r r dr Rt t    t              00 . . . ln ln . . .. r l c l l f R r d R r r dr Rt t t    t                2 2 . ln . ln . .ln . . 2 . 2 r l l R r R r A R R r rdr R R   t                         (5) Với   2 0 . . , . . l c l f r A s m t t         Đặt .ln . R R I r r dr  (6) Tích phân từng phần (6) với lnu r và .dv r dr có kết quả    2 2 2 21 1. .ln .ln . 2 4 I r r R R r R    . Thay kết quả (6) vào (5) ta có:      2 2 2 21. ln ln ln 2 . 2 lAr R r R r r R R R  t                (7) Khi kết thúc quá trình chuyển pha  2 0rt t      2 2 0 . . . . 1 2. , .4. . l c l l f r R s Rt t    t           (8) 2.3.3. Tính thời gian quá lạnh băng 3t theo trình tự sau a. Mô tả trường nhiệt độ ( , )t r t như Hình 2 và quá trình quá lạnh 3t (, t) = ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎧ 1.ℎ t = 0, (, t = 0) = t = , ∀ ∈ (0, ) 2. ℎ t ∈ (0, t), (, t) = − − . .. .. 3. ℎ t = t = . . 2. . − − ℎì (, t) = , ∀ 4. ℎ t ∈ (t, t + t)ℎì (, t) = , ∀ ∈ (0, ()) (t) + (t), ∀ ∈ ( ÷ ) 5. ℎ t = t + t = t ℎì (, t) = ạ = 0 () + (), () = () 6. ℎ t ∈ (t + t = 0, (t + t + t) ℎì (, t) = (t) + (t) 7. ℎ t = = , ℎì (, t) = () + () b. Lập công thức tính 3t Phương trình cân bằng nhiệt cho khối băng 2. .V R l trong khoảng thời gian quá lạnh 30t t     là: ISSN 1859-1531 - TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ ĐẠI HỌC ĐÀ NẴNG, SỐ 11(96).2015, QUYỂN 2 13       3 30 . . .2 . .fU U U t t r lt t  t     w (9) Suy ra         33 0 , . .f U U s t t F t t  t   w , với các thông số      30 , ,U U tt t  w xác định theo phân bố        2 2 3 3, , 0, , , 0t r r Rt t t t t t      phần tiếp theo. c. Xác định gần đúng phân bố  ,t r t khi quá lạnh: Trong quá trình quá lạnh,    3 1 2 1 2 30 ,t t t t t t t     thì phân bố chính xác của  ,t r t ,      1 3 2 3, . .ln .t r C r Ct t t  đi qua hai điểm  0 0, 0,t t  và / ,s fR R t     với    0 . 0, .exp 2. . . f f k t t t t cR t t      . Vì hàm ln rkhông khả vi tại r = 0 nên để tính gần đúng, có thể coi đoạn cong 1 2 1 2.ln .C r C C r C   , tức là tuyến tính hóa đường cong  ,t r t . Khi đó phân bố  ,t r t trong quá trình quá lạnh có dạng:      1 2, .lnt r C r Ct t t  đi qua hai điểm  0 0, (0, )r t t t  và / ,s fR r R t t       với  1C t và  2C t xác định theo:                       0 1 2 1 0 1 2 2 0 .0 / / . / f s f s t C C C t t R t C R C C t t t t t t   t   t t t                   (10) Trường  ,t r t khi quá lạnh, gần đúng là:      0 0, . / f s t t t r r t R t t t       (11) d. Tính số gia nội năng của V trong quá trình quá lạnh 3t Nội năng tức thời lúc t của 2. . . .dV r l dr tại r khi 1 2.t C r C  là:    1 2. . . , 2. . . . . . . .s s s sdU dV c t r r l c dr C r C t     Tổng nội năng U của 2. .V R l lúc t là:      3 21 2 1 2 0 2. . . . . . . 2. . . . . . . /3 . /2 R s s s s r U l c C r C rdr rl c C R C R    t t               2 1 2 . . . . . 2. 3. 3 s sc l RU C C   t t t     Số gia U của V trong quá trình quá lạnh 3t là:    30U U U t             2 1 1 3 2 2 3 . . . . . 2 . 0 3. 0 3 s sc l RU R C C C C   t t      (12) Thay các trị số  1 0C ,  1 3C t ,  2 0C ,  2 3C t từ (9) vào (11), sẽ có:       2 0. . . . . 3. . . , 3 s s s cc l R R t t U J        (13) e. Tính nhiệt độ trung bình của mặt trụ  t tw trong thời gian  30t t   : Nhiệt độ mặt trụ lúc t khi quá lạnh là  t tw xác định theo quan hệ:           0 0 / / 1 . / / f s f f s f s t t R t t t t R t t t   t t   t            w w Lấy trung bình theo  30t t   có    0 0 / 2ct t tt   và      0 2 / . / 2f c f s st t t t t Rt         w (14) f. Xác định công thức tính thời gian quá lạnh 3t Theo (8) ta có:        03 0 . . . 3 . , 3. . . 2 s s s c s c f c R R t t s t t t    t         (15) 2.4. Tính toán thời gian cấp đông vật ẩm dạng cầu 2.4.1. Tính thời gian làm lạnh vật ẩm 1t Thiết lập tương tự như mục 2.3.1 cho vật ẩm dạng cầu có thể tích 3 4 . . 3 cau cauV R ta có biểu thức:  1 0 . . .ln , . i fcau l l cau f t tV c s F t t  t      (16) 2.4.2. Tính thời gian chuyển pha 2t Gọi r là bán kính lớp băng đã tạo ra trước thời điểm t , dr là lớp băng mới tạo ra sau thời gian dt . Phương trình cân bằng nhiệt tức thời lúc t cho khối cầu băng mới tạo ra dV = Fcau.dr = 4.π.r2.dr sau thời gian vô cùng bé dt là:  02 2 4. . . . . . 4 . 1 1 1 1 . . f l c l cau caucau t t r dr r d r RR     t            (17)   2 2 0 . . 1 1 1 1 . . . . . l c cau l cau caucauf r d r dr r RRt t   t                   0 2 2 00 . . 1 1 1 1 . . . . . cau cau r l c cau l cau caucaufR r d r dr r RRt t t   t                       2 2 0 . . . . 1 2. , .6. . l c cau l caul f r R s Rt t    t           (18) 2.4.3. Tính thời gian quá lạnh băng 3t Thiết lập tương tự như mục 2.3.3 cho vật ẩm dạng cầu có thể tích 3 4 . . 3 cau cauV R và diện tích 24. . caucau F R , ta có công thức tính thời gian quá lạnh 3t :        03 0 . . . 4 . , 6. . . 2 s s cau s cau c s c f c R R t t s t t t    t         (19) 2.5. Khảo sát bài toán cấp đông thực phẩm và so sánh với kết quả thực nghiệm 2.5.1. Xác định hệ số tỏa nhiệt bề mặt  Hệ số tỏa nhiệt bề mặt giữa thực phẩm và không khí làm lạnh đối với vật ẩm dạng trụ không bao gói có thể xác định bằng công thức [4]: 0,6 212,5. , / .tbv m K    W (20) Để đảm bảo độ chính xác của nghiên cứu thực nghiệm, tác giả đã tiến hành đo tốc độ gió của buồng cấp đông IQF 14 Nguyễn Bốn, Võ Chí Chính, Hoàng Minh Tuấn tại các vị trí khác nhau. Kết quả tốc độ gió trung bình được sử dụng trong tính toán 2,8 /tbv m s . 2.5.2. Khảo sát bài toán cấp đông thực phẩm Cần cấp đông cho khối xúc xích dạng hình trụ có 0,015D m và 0,15l m bằng tủ đông IQF thẳng. Nhiệt độ ban đầu của xúc xích ti = 31oC được làm lạnh đến nhiệt độ tâm theo yêu cầu tc = -9,5oC. Nhiệt độ của không khí lạnh là tf = -38oC. Tính toán thời gian cấp đông. Bảng 1. Thành phần khối lượng của xúc xích heo hầm hạt sen Bé Khỏe Công ty Thực phẩm Gia đình Anco STT Thành phần Ký hiệu Giá trị và đơn vị 1 Nước W 62,9% 2 Béo F 10,45% 3 Rắn S 26,65% (Phân tích mẫu xúc xích tại TT kiểm nghiệm thuốc Mỹ phẩm và Thực phẩm, Sở Y Tế tỉnh Thừa Thiên Huế) Bảng 2. Các thông số vật lý của xúc xích heo hầm hạt sen Bé Khỏe được tính toán theo thành phần khối lượng STT Thông số Ký hiệu Giá trị và đơn vị Chú thích Pha ẩm (l) Chuyển pha (l) Pha rắn (s) 1 Khối lượng riêng  1045kg/m3 1045kg/m3 1045kg/m3 [4] 2 Nhiệt độ bắt đầu đóng băng 0 t -1,171 oC [4] 3 Nhiệt ẩn hóa rắn cr 228kJ/kg [2] 4 Nhiệt dung riêng c 3200,9J/kg.K 2014,9J/kg.K [4] 5 Hệ số dẫn nhiệt  0,4634W/m.K 1,1195W/m.K [4] Bảng 3. Kết quả tính toán thời gian đóng băng Thời gian Công thức tính cho vật dạng trụ Trị số, [s] Công thức tính cho vật dạng cầu Trị số, [s] 1t 1 0 . . .ln . i fl l f t tV c F t t  t     340 1 0 . . .ln . i fcau l l cau f t tV c F t t  t     227 2t   2 2 0 . . . . 1 2. .4. . l c l l f r R Rt t    t          782   2 2 0 . . . . 1 2. .6. . l c cau l caul f r R Rt t    t          521 3t       0 3 0 . . . 3 . 3. . . 2 s s s c s c f c R R t t t t t    t        91       0 3 0 . . . 4 . 6. . . 2 s s cau s cau c s c f c R R t t t t t    t        60 ltt 1 2 3ltt t t t   1213 1 2 3ltt t t t   808 2.5.3. Thực nghiệm kiểm chứng Để so sánh kết quả lý thuyết với thực nghiệm tác giả đã tiến hành. a. Chọn vật ẩm, thiết bị đo và hệ thống lạnh + Vật ẩm là xúc xích heo hầm hạt sen Bé Khỏe có các thông số vật lý như Bảng 2. + Thiết bị đo gồm: đồng hồ đo thời gian có độ chính xác 0,01s, nhiệt kế thermo scientific với độ chính xác 0,1oC và máy đo tốc độ gió Prova AVM 03 sai số 3% . + Hệ thống lạnh: IQF thẳng tại Công ty cổ phần Thủy sản Sông Hương - Thừa Thiên Huế. b. So sánh thời gian cấp đông thực nghiệm và thời gian lý thuyết cho vật ẩm dạng trụ với cách tính sai số  t theo công thức    / , %tn lt tnt t t t  Bảng 4. Kết quả thực nghiệm thời gian cấp đông TT , o it C   , oft C   , oct C    ,lt pt  ,tn pt  , %t 1 31 -38 -9,5 20,22 21,90 7,67 2 31 -38 -12,5 20,86 23,10 9,67 3 12 -38 -18 19,32 21,30 9,29 4 31,2 -36 -12,8 22,10 24,20 8,68 5 18 -36 -13,3 20,25 22,10 8,37 Sai số trung bình thực nghiệm 8,74 c. Một số hình ảnh thực nghiệm Hình 3. Đo nhiệt độ ban đầu và tâm của thực phẩm Hình 4. Thực phẩm đi vào và ra hệ thống IQF ISSN 1859-1531 - TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ ĐẠI HỌC ĐÀ NẴNG, SỐ 11(96).2015, QUYỂN 2 15 Hình 5. Thiết bị đo tốc độ gió trong tủ IQF 3. Kết luận và kiến nghị - Đây là phương pháp đơn giản để dự đoán thời gian đóng băng cho vật ẩm có dàng hình trụ vô hạn và cầu. Mô hình này cho phép tính toán đơn giản và nhanh chóng dựa trên phương trình cân bằng nhiệt; - So sánh ( 1 2 3, ,t t t ) của vật ẩm dạng trụ vô hạn và cầu thì thời gian làm lạnh, chuyển pha và quá lạnh cho vật ẩm dạng trụ vô hạn nhiều hơn khoảng 1,5 lần; - Sai số giữa lý thuyết và thực nghiệm trung bình không quá 10%, cho thấy: các công thức đưa ra trong phương pháp này khá chính xác và có thể ứng dụng để dự đoán quá trình đóng băng thực phẩm trong thực tế; - Các sai số gây ra chủ yếu trong phương pháp này do: các thông số vật lý ( , , ...c  ) chưa có số liệu chính xác, ảnh hưởng của thiết bị và điều kiện môi trường thực nghiệm. TÀI LIỆU THAM KHẢO [1] Nguyễn Bốn, Tính toán thời gian đông lạnh thực phẩm, Tạp chí Khoa học và Công nghệ Nhiệt, số 7/2004. [2] Becker, B.R. and B.A. Fricke. (1999). Food thermophysical property models. International Communications in Heat & Mass Transfer 26(5):627-636. [3] Castaigne, F. (1985b). Calcul des temps de congelation d'aliments ayant la forme d'un cylindre infini, d'un cylindre fini ou d'une sphere. Lebensm. Wiss. u. Technol. 18:137. [4] Cleland, DJ., Cleland, AC., White, SD., Love, RJ., Merts, I., East, A., Paterson, AHJ.(2010). Cost-Effective Refrigeration. Palmerston North, New Zealand: Massey University. [5] Cleland, A.C. and Earle, R.L. (1977a). A comparison of analytical and numerical methods for predicting the freezing times of foods. 1. Food Sci. 42:1390. [6] Cleland, A.C. and Earle, R.L. (1977b). The third kind of boundary condition in numerical freezing calculations. Int. J. Heat Mass Transfer. 20: 1029. [7] Cleland, A.C. and Earle, R.L. (1982). A simple method for prediction of heating and cooling rates in solids of various shapes. Rev. Int. du Froid, 5:98. [8] Cleland, A.C. and Earle, R.L. (1987a). Prediction of freezing & thawing times for multidimensional shapes by simple formulae. Part 1: Regular shapes. International Journal of Refrigeration, 10, 157-164. [9] Cleland, A.C. and Earle, R.L. (1987b). Prediction of freezing & thawing times for multidimensional shapes by simple formulae. Part 2: Irregular shapes. International Journal of Refrigeration, 10, 234-240. [10] De Michelis, A. and Calvelo, A. (1982). Mathematical models for symetric freezing of beef. J. Food Sci, 47:1211. [11] De Michelis, A. and Calvelo, A, (1983). Freezing time predictions for brick and cylindrical shaped foods. J. Food Sci. 48:909. [12] Hossain, Md.M., Cleland, D.J., Cleland, A.C. (1992b). Prediction of freezing and thawing times for foods of two-dimensional irregular shape by using a semi-analytical geometric factor. International Journal of Refrigeration, 15, 235-240. [13] Hung, Y.C. and Thompson, D.R. (1983). Freezing time prediction for slab shape foodstuffs by an improved analytical method. J. Food Sci. 48:555. [14] McNabb, A., Wake, G.C., Hossain, Md.M., Lambourne, R.D. (1990b) Transition times between steady states for heat conduction, Part II: Approximate solutions and examples. Occasional Pubs in Maths & Statistics No.21, Massey University. [15] Nagaoka, J., Takaji, S., and Hohani, S. (1955). Experiments on the Freezing of fish in air blast freezer, Proc. IX Int. Congo Refrig. 4: 105. (BBT nhận bài: 17/08/2015, phản biện xong: 10/09/2015)