BỘ GIÁO DỤC VÀ ĐÀO TẠO
TRƯỜNG ĐẠI HỌC BÁCH KHOA HÀ NỘI

ĐỖ HỮU HOÀNG

NGHIÊN CỨU MÔ PHỎNG VÀ XÁC ĐỊNH CHẾ ĐỘ CẤP
ĐÔNG HỢP LÝ CHO CÁ TRA VIỆT NAM
Chuyên ngành:: Kỹ Thuật Nhiệt
Mã số:

62520115

LUẬN ÁN TIẾN SĨ KỸ THUẬT NHIỆT

NGƯỜI HƯỚNG DẪN KHOA HỌC:
ĐẶNG QUỐC PHÚ

1. GS.TSKH.

:

2. TS.

: NGUYỄN VIỆT DŨNG

Hà Nội – 2014


CHƯƠNG 1
TỔNG QUAN VỀ LĨNH VỰC NGHIÊN CỨU
1.1

triển nhanh như sản phẩm cá tra (và basa) của Việt Nam xem hình 1.1 và hình 1.2. Trong
vòng 10 năm qua, sản lượng cá tracủa Việt Nam đã tăng 50 lần, giá trị xuất khẩu tăng 65%
lần và hiện đang chiếm tới 90% thị phần thế giới [1-3,15].

Cá đông lạnh

Giáp xác

Thân mềm

300
250
200
150
100
50
0
Năm

Hình 1.1 Kim ngạch xuất khẩu cá đông lạnh và giáp xác, thân mềm của Việt nam 2002 2011 (triệu USD)[1-3]
1


Hình 1.2 Tình hình xuất khẩu cá tra từ năm 2009 đến ngày 15/6/2013[1-3]
Ở nước ta khu vực sản xuất cá da trơn là đồng bằng sông Cửu Long xuất khẩu sản
phẩm cá tra đem lại nguồn thu to lớn và quan trọng cho đất nước bởi vậy duy trì và phát
triển bền vững nguồn lợi từ xuất khẩu cá da trơn là rất quan trọng.

1.1.2 Quy trình công nghệ chế biến và bảo quản đông lạnh cá da trơn tại khu
vực đồng bằng sông Cửu Long - Việt nam

Vì lý do nêu trên, trong khuôn khổ luận án này, đối tượng nghiên cứu là quá trình
cấp đông cá tra bằng các thiết bị cấp đông rời dạng IQF.
2


Hình 1.3 Quy trình công nghệ chế biến và bảo quản tôm và cá fillet[15]

1.1.3 Tiêu hao năng lượng trong chế biến cá da trơn
Điện năng là năng lượng tiêu thụ chính trong các nhà máy chế biến cá da trơn tại
khu vực đồng bằng sông cửu long, theo thống kê tại các nhà máy chế biến, suất tiêu thụ
điện năng để chế biến 1kg các da trơn từ nguyên liệu đến thành phẩm dao động từ 0.4 đến
0,45 kWh/kgSP. Số liệu này tính trung bình chung, tuy nhiên giá trị này rất biến động tùy
thuộc và sản lượng nguyên liệu cung cấp cho nhà máy. Cụ thể tiêu hao điện năng điển hình
của thiết bị trong nhà máy chế biến được trình bày trên hình 1.4 [15]
Tỷ lệ của các hộ theo tiêu thụ điên
Chiếu sáng
2%

Văn phòng
2%

Thiết bị
sản xuất
96%
Máy nén
72.86%

Tỷ lệ theo tiêu thụ điện
Thiết bị phụ trợ
6%

Trạm bơm cấp
nước+ xử lý nước
5.66%

Quạt dàn lạnh
2.50%

Hình 1.4 Tỷ lệ tiêu thụ điện năng điển hình của các thiết bị tiêu thụ điện tại các nhà máy
[15]
Theo thống kê về suất tiêu hao năng lượng [18],[50],[51],[60],[61],[97] đánh giá
suất tiêu hao năng lượng trong chế biến thủy sản giữa Việt nam, các nước đang phát triển
và các nước phát triển như sau:
3


Bảng 1.1 Tiêu hao năng lượng trong chế biến thủy sản
Stt Ngành
1

Suất tiêu
hao điện

Loại SP

Thủy
sản

637,11

Đơn vị

Tỷ lệ chênh lệch
suất tiêu hao năng
lượng giữa VN
và thế giới(%)

1.166,4-15.883,2

68

Từ các số liệu ở trên chúng ta thấy rằng để đảo bảo khả năng cạnh tranh sản phẩm
thủy hải sản đông lạnh nói chung và cá da trơn nói riêng việc giảm tiêu hao năng lượng
trong quá trình chế biến và bảo quản sản phẩm thủy hải sản đông lạnh là vấn đề cấp thiết,
đặc biệt là trong giai đoạn làm lạnh và cấp đông trong giai đoạn này năng lượng sử dụng
chiếm trên 70%.

1.1.4 Hao hụt sản phẩm trong quá trình cấp đông
Đối với sản phẩm thủy hải sản trong quá trình đông lạnh, đặc biệt là các sản phẩm
có giá trị kinh tế cao hao hụt khối lượng trong quá trình cấp đông ảnh hưởng rất lớn đến
giá thành sản phẩm. Hao hụt khối lượng sản phẩm trong cấp đông do bốc hơi nước từ bề
mặt sản phẩm vào không khí do có sự chênh lệch của áp suất riêng phần của hơi nước
trong quá trình cấp đông.
Đối với công nghệ cấp đông cá da trơn bằng IQF độ hao hụt khối lượng sản phẩm
tỷ lệ với diện tích bề mặt sản phẩm và thời gian cấp đông. Khi tiết diện bề mặt sản phẩm
càng lớn thì độ hao hụt càng lớn. Ảnh hưởng của thời gian cấp đông đến độ hao hụt khối
lượng sản phẩm cấp đông được thể hiện trên đồ thị sau:

1,0

0,5


tiêu chuẩn, do đó việc xác định chính xác tâm của miếng cá là không khả thi. Đồng thời
định nghĩa tâm cũng không rõ ràng là tâm thấm nhiệt hay là tâm hình học. Do đó kết quả
xác định thời gian cấp đông là hoàn toàn chủ quan phụ thuộc vào người thực hiện và kinh
nghiệm của họ. Việc này ảnh hưởng không nhỏ tới chất lượng sản phẩm và tiêu hao năng
lượng cho quá trình cấp đông sản phẩm. Do đó cần thiết phải có những nghiên cứu để dự
đoán thời gian cấp đông của sản phẩm chính xác hơn.

1.1.6 Các yếu tố ảnh hưởng tới chất lượng cấp đông cá
Đối với sản phẩm thủy hải sản nói chung và cá da trơn nói riêng yếu tố quyết định
chất lượng sản phẩm là tốc độ cấp đông. Tốc độ càng nhanh, thời gian cấp đông càng ngắn
thì chất lượng sản phẩm càng tốt và thời gian bảo quản càng được kéo dài. Tuy nhiên để
thực hiện được điều này, nhiệt độ môi trường cấp đông phải rất thấp, tốc độ gió cao, dẫn
tới công suất điện tiêu thụ của hệ thống lạnh tăng lên, hiệu suất năng lượng giảm đi, giá
thành sản phẩm tăng cao. Do đó, việc đánh giá ảnh hưởng của các thông số môi trường cấp
đông là nhiệt độ và tốc độ gió có ý nghĩa rất quan trọng trong việc kiểm soát chất lượng
sản phẩm và tiêu hao năng lượng trong quá trình làm lạnh và cấp đông.
Trên thực tế sản xuất nghiên cứu vấn đề này bằng thực nghiệm là hết sức khó khăn
và thậm chí là bất khả thi, do chiều dài buồng cấp đông lớn và có cấu trúc đóng kín với
trường nhiệt độ và tốc độ là thông số rải, không cho phép chúng ta đo đạc trực tiếp khi
thiết bị đang vận hành. Hơn nữa cần thiết phải thực hiện một khối lượng thí nghiệm không
hề nhỏ để có thể giải quyết vấn đề trên, đòi hỏi chi phí rất lớn. Do vậy giải pháp tốt nhất là
xây dựng mô hình mô phỏng quá trình làm lạnh và cấp đông cá da trơn để nghiên cứu ảnh
hưởng của các thông của môi trường cấp đông. Mô hình này cho phép chúng ta xác định
các thông số (1) thời gian cấp đông, (2) trường nhiệt độ của thực phẩm.

1.1.7 Kết luận
Cá da trơn là nguồn lợi thủy sản lớn của Việt Nam với tổng kim ngạch xuất khẩu
đạt gần hai tỉ USD trong thời gian gần đây. Trong đó cá tra chiếm chủ đạo tới khoảng 90%
Tuy nhiên việc phát triển xuất khẩu mặt hàng cá da trơn đông lạnh đang gặp nhiều thách
thức. Trong đó thách thức lớn là chất lượng chế biến chưa cao do công nghệ cấp đông chưa

1.2 Tổng quan về mô phỏng quá trình làm lạnh cấp đông thực phẩm
1.2.1 Vai trò của việc nghiên cứu mô phỏng quá trình lạnh đông
Trong quá trình làm lạnh hay cấp đông rất quan trọng: (1) dự đoán được thời gian
cấp đông, (2) trường nhiệt độ của thực phẩm theo thời gian. Giải quyết được vấn đề (1)
giúp giải quyết bài toán nâng cao chất lượng chế biến, giải quyết được vấn đề (2) cho
phép thiết kế, chọn lựa hệ thống lạnh hợp lý (bài toán tiết kiệm năng lượng). Muốn dự
đoán được các thông số trên và đánh giá chất lượng quá trình chế biến lạnh trong điều kiện
thực nghiệm không phải bao giờ cũng thực hiện được, cùng với với sự đa dạng về đối
tượng cấp đông và thiết bị trong thực tế sản xuất. Rất cần thiết phải kết hợp lý thuyết - thực
nghiệm xây dựng được mô hình làm lạnh, cấp đông cho thực phẩm. Đây là vấn đề được
quan tâm cả ở trong và ngoài nước từ mấy chục năm trở lại đây.

1.2.2 Thực trạng của việc nghiên cứu mô phỏng quá trình lạnh đông
Cho tới nay đã có hàng chục dạng mô hình dùng để mô phỏng quá trình cấp đông
thực phẩm được đưa ra sử dụng trong và ngoài nước. Tuy nhiên về bản chất các mô hình
này đều dựa trên cơ sở giải hệ phương trình vi phân dẫn nhiệt phi tuyến kết hợp với điều
kiện biên cho thực phẩm ở môi trường cấp đông. Phương trình vi phân dẫn nhiệt tổng quát
viết cho một phân tố của vật thể được cấp đông có dạng như sau:
∂T(r, τ)
C(T )ρ(T )
= div[λ(T ) ⋅ gradT(r, τ)] + q v (r, τ)
(1.1a)
∂τ
Trong đó:
C(T)- nhiệt dung riêng phụ thuộc vào nhiệt độ của thực phẩm, kJ/kg.K
ρ(T)- khối lượng riêng phụ thuộc vào nhiệt độ của thực phẩm, kg/m3
λ(T)- hệ số dẫn nhiệt phụ thuộc vào nhiệt độ của thực phẩm, W/m.K
qv(r, τ)- nguồn nhiệt trong, sinh ra do sự chuyển pha của nước phụ thuộc vào tọa độ
của phân tố và thời gian τ, W/m3
T(r, τ)- nhiệt độ của phân tố phụ thuộc vào tọa độ và thời gian, K.


r
r
T(x, τ = 0) = Tin (x) ,

(1.1c)
Lời giải của hệ (1.1) cho chúng ta phân bố trường nhiệt độ trong lòng sản phẩm và
thời gian cấp đông.
Tuy nhiên khi giải hệ phương trình vi phân dạng (1.1) gặp phải một số vấn đề khó
khăn. Tính chất nhiệt vật lý của thực phẩm như nhiệt dung riêng, hệ số dẫn nhiệt thay đổi
đột ngột lân cận điểm đóng băng, dẫn đến những phương trình vi phân từng phần có độ phi
tuyến cao, rất phức tạp để giải. Đối với những vật thể có hình dạng phức tạp như thực
phẩm, quá trình đóng băng càng khó dự đoán. Lạnh đông trong thực tế bao gồm một vài
hiện tượng vật lý diễn ra đồng thời: truyền nhiệt, truyền chất, sự lớn lên của mầm tinh thể,
thay đổi thể tích, căng cơ học và các ứng suất.
Các cách tiếp cận khác nhau để giải quyết hệ phương trình vi phân dạng (1.1) sẽ
cho các dạng mô hình mô phỏng khác nhau của quá trình cấp đông thực phẩm với độ chính
xác khác nhau. Như đã trình bày ở trên các dạng mô hình này rất đa dạng, tuy nhiên nhìn
chung có thể chia làm hai dạng chính liên quan tới phương pháp giải quyết bài toán nêu
trên. (1) các lời giải dựa trên phương pháp giải tích, (2) các lời giải dựa trên phương pháp
số. Dưới đây chúng ta sẽ khảo sát chi tiết từng phương pháp.

1.2.3 Phương pháp giải tích
Phương pháp giải tích dựa trên giả thuyết quá trình chuyển pha là lý tưởng, sự
chuyển pha và giải phóng nhiệt ẩn đóng băng diễn ra ở nhiệt độ điểm băng không đổi,
thường được ký hiệu là Tf, đồng thời tồn tại bề mặt phân pha giữa vùng đóng băng và vùng
chưa đóng băng. Bỏ qua thành phần nhiệt hiện giải phóng trong quá trình kết đông, thông
số nhiệt vật lý được xem là hằng số, thực phẩm được xem là dung dịch đồng chất và đẳng
hướng. Kết quả nhận được bằng cách giải mô hình 1.1a trong hai vùng tương ứng với cùng
điều kiện biên T = Tf ở bề mặt phân pha.

vật có hình dáng xác định nhiều chiều (thanh hình chữ nhật dài vô hạn, hình trụ hữu
hạn…), biểu thức giải tích cho hệ số E cho dưới dạng chuỗi vô hạn (McNabb, 1990a,
1990b)[69-70] và đồ thị (Hossain, 1992a) [71].Nhìn chung các hệ số này được xác định
bằng phương pháp hồi quy từ số liệu thực nghiệm cho trước [19-22, 62-63,106-108..].
Bảng 1.3 trình bày tổng hợp các dạng công thức tính toán thời gian cấp đông là biến thể
của dạng nghiệm của phương trình Plank (1.2).

Bảng 1.3 Tổng hợp các mô hình xác định thời gian cấp đông điển hình
Stt

Tác giả

1

R Plank

2

J Nagaoka, S Takagi,
S Hotan

3

FL Levy

Mô hình toán
τ Plank =

ρL f  2PR 4QR 2 
+



λi E
 N Bi Nste Nste 
1,65Nste  t − t f
0 , 0017N Bi + 0 ,1727N pk 
×10
ln
1 −
λi
 t ref − t f


τ slab =
5

Q.T. Phạm





τtam phang
ρa  ∆H 1 ∆H 2   N Bi  τ
=
+
1
+
;
batky

8


Từ bảng 1.3 có thể thấy những hạn chế của phương pháp giải tích là chỉ cho thời
gian cấp đông, mà không xác định được trường nhiệt độ.
Hơn nữa những dạng công thức trong bảng 1.3 chỉ có thể áp dụng cho một số loại
thực phẩm có thông số nhiệt vật lý được xác định trước. Đồng thời cũng một câu hỏi đặt ra
là các đại lượng nhiệt vật lý được sử dụng trong các công thức ở trên như hệ số trao đổi
nhiệt đối lưu α, hệ số dẫn nhiệt λ, hiệu enthalpy của thời điểm đầu và cuối của quá trình
cấp đông ∆H sẽ được tính toán như thế nào nếu không biết được trường nhiệt độ, hàm
lượng nước đóng băng phụ thuộc vào nhiệt độ.
Giả thiết quá trình chuyển pha là lý tưởng, sự chuyển pha và giải phóng nhiệt ẩn
đóng băng diễn ra ở nhiệt độ điểm băng không đổi đồng thời tồn tại bề mặt phân pha giữa
vùng đóng băng và vùng chưa đóng băng là không đúng. Vì trong khi thực tế trong quá
trình kết đông diễn ra ở nhiệt độ thay đổi do sự phát triển của mầm tinh thể băng và không
tồn tại bề mặt phân pha
Do đó với cá tra tại Việt nam đặc biệt với thông số nhiệt vật lý chưa được xác định,
không thể áp dụng các mô hình toán trên để xác định thời gian cấp đông của thực phẩm.
Giải pháp tốt nhất để mô phỏng quá trình cấp đông cá tra là giải phương trình dạng (1.1)
bằng phương pháp số.

1.2.4 Phương pháp số giải hệ phương trình vi phân dẫn nhiệt
Bản chất phương pháp số để giải hệ phương trình (1.1) bao gồm hai bước: rời rạc
hóa các miền liên tục để thu được một bộ phương trình vi phân thường (ODE) tương ứng
với các nút nhiệt độ, sau đó giải bộ phương trình ODE bằng cách chuyển phương trình vi
phân về phương trình đại số tuyến tính qua phép xấp xỉ sai phân. Bộ phương trình vi phân
thường có thể viết dưới dạng ma trận như sau:

C



Tip, j+1 − Tip, j
∆t

=

λ  Ti −1, j − 2Ti , j + Ti +1, j Ti , j−1 − 2Ti , j + Ti , j+1 
+

C.ρ 
∆x 2
∆y 2


(1.4)

Nhiệt độ các điểm nút (i,j) được xác định từ lời giải 2 sơ đồ sau:
- Sơ đồ hiện (sai phân tiến)
p +1
i, j

T

 Tip−1, j − 2Tip, j + Tip+1, j Tip, j−1 − 2Tip, j + Tip, j+1 
λ
p
=
∆t 
+
 + Ti , j

Phương pháp thể tích hữu hạn tập trung vào điểm giữa phân tố thể tíchtương tự như
phương pháp SPHH (Malanvà cộng sự 2002).
Nhiệt độ tại các điểm nút được xác định tương tự phương pháp sai phân hữu hạn, bằng
cách rời rạc phương trình ma trận đặc trưng theo thời gian ứng với thời gian p (sai phân
lùi) hoặc p+1 (sai phân tiến), cụ thể như sau. Áp dụng định luật bảo toàn năng lượng, tại tất
cả các phân tố cũng như trong toàn vật thể ta có:



∂T

∫ ρc ∂τ

Ω


-∇ ( λ∇T ) - q dΩ = 0


(1.7)

Sử dụng biến đổi tích phân thể tích về tích phân mặt phương trình 1.7 trở thành:
 ∂T

(1.8)
∫Ω ρc ∂τ - q dΩ - ∫S n ( λ∇T )dS = 0
Với S là diện tích bề mặt truyền nhiệt của phân tố, n là phương pháp tuyến bề mặt
truyền nhiệt và Ω là phân tố thể tích. Thành phần đầu của phương trình 1.8 là độ tăng
enthalpy của phân tố có thể tính xấp xỉ theo công thức:


Với Bij là hệ số phụ thuộc vào cách sắp xếp các phân tố thể tích.
Tổng hợp các phương trình tại điểm nút ta sẽ thu được phương trình ma trận có dạng 1.11.

δV.ρ.c

∂Ti n
= ∑ λBijTj + qδV
∂t i =1

(1.11)

- Phương pháp Phần tử hữu hạn (PTHH) là phương pháp số để giải các bài toán
được mô tả bởi các phương trình vi phân đạo hàm riêng cùng với các điều kiện biên cụ thể.
10


Cơ sở của phương pháp này là rời rạc miền nghiệm liên tục và phức tạp của bài toán thành
các miền con gọi là các phần tử hữu hạn. Tuỳ theo yêu cầu của bài toán mà các miền con
tức các phần tử hữu hạn này có cấu trúc khác nhau, tinh xảo và liên kết với nhau bởi các
nút. Việc tìm lời giải chính xác của bài toán được thay thế bằng việc tìm dạng gần đúng tại
các nút thông qua hàm xấp xỉ trên từng phần tử. Hàm xấp xỉ được gọi là phương trình đặc
trưng của phần tử, có thể được xác định bằng phương pháp biến phân, hoặc phương pháp
số dư trọng số.
Theo phương pháp PTHH, nhiệt độ tại các điểm nút được xác định như sau:
- Từ phương trình vi phân dẫn nhiệt, chuyển về phương trình ma trận đặc trưng:

[C] ∂T  + [K ]{T} = {f }

 ∂τ 
T

 1
[C] + [K ]{T p+1} = {f } + [C] 1 {T p }+ 1 − θ {T p }

1
 θ∆τ

 θ∆τ


(1.13)

Khi θ =0 (sơ đồ hiện hoàn toàn):

Tp+1 = [C] ∗ [([C] − [K]∆τ){TP } + {f P }∆τ]
−1

Khi θ =0 (sơ đồ ẩn hoàn toàn):
−1
T p+1 = [[C] − [K ]∆τ] ∗ [[C]T p + {f P +1}∆τ]

(1.14)
(1.15)

Khi θ =1/2, sơ đồ nửa ẩn

{[C] + 0,5∆τ[K]}{T}p+1 ={[C] − 0,5∆τ[K]}{T}p + 0,5∆τ({f}p+1 + {f}p )

(1.16)
Tóm lại phương pháp số về nguyên tắc có thể cho lời giải với độ chính xác yêu cầu
trong trường hợp tổng quát, mặc dù vậy trong thực tế độ chính xác của chúng bị giới hạn

(1) mô hình mô phỏng quá trình lạnh đông bằng phương pháp giải tích và công thức thực
nghiệm chỉ yêu cầu mô hình nhiệt vật lý đơn giản hoặc không đầy đủ (phương trình xác
định thời gian lạnh đông của Plank, [92], các công thức thực nghiệm của Phạm [82-90],
Clealand & Earle[19-24], Salvadori & Mascheroni [93]), (2) một số mô hình phương pháp
số không có khả năng giải cho những mô hình nhiệt vật lý phức tạp, tính chất nhiệt vật lý
là hàm của nhiệt độ và thành phần thực phẩm – bài toán dẫn nhiệt phi tuyến (mô hình của
Zengfu Wang, Han Wu, Guanghua Zhao, Xiaojun Liao, Fang Chen, Jihong Wu, Xiaosong
Hu, 2007).
Tocci, A.M & Mascheroni, R.H (1995)[79] trong một nghiên cứu về mô hình xác
đinh thời gian lạnh đông thịt bò băm viên sử dụng 3 mô hình nhiệt vật lý khác nhau của
Sanz (1989), Cleland & Earle (1984) và Cleland & Earle (1986) đã chỉ ra rằng: (a) dữ liệu
về tính chất nhiệt vật lý phù hợp là chìa khóa đạt được độ chính xác trong mô hình dự đoán
thời gian lạnh đông, (b) tính chất nhiệt vật lý phải kể đến ảnh hưởng của nhiệt độ và thành
phần thực phẩm, (c) mỗi một mô hình mô phỏng được phát triển dựa trên một mô hình
nhiệt vật lý nhất định và kết quả dự đoán sẽ chính xác hơn nếu mô hình đó được sử dụng.
Trong các nghiên cứu về mô phỏng tính chất nhiệt vật lý của thực phẩm cách tiếp cận của
trường phái Latushev và Chumark I.G & Onistchenko V.P.[35-37], dựa trên cơ sở coi thực
phẩm như một hệ nhiệt động đa thành phần, cân bằng và có chuyển pha cho kết quả khả
quan và có cơ sở lý thuyết chặt chẽ hơn cả. Dưới đây chúng ta sẽ lần lượt xem xét các mô
hình nêu trên.

1.3.1

Thành phần băng

Trong thực phẩm, nước chiếm tỷ trọng rất lớn. Sự chuyển pha của thành phần nước từ
lỏng thành băng quyết định rất lớn đến sự thay đổi tính chất nhiệt vật lý của thực phẩm vì
tính chất của nước và băng là hoàn toàn khác nhau. Dự đoán chính xác thành phần băng
trong thực phẩm theo nhiệt độ đóng vai trò quan trọng trong việc xác định những đặc điểm
nhiệt vật lý và sự biến đổi enthanlpy của thực phẩm trong quá trình lạnh đông, bảo quản và


Phân tử lượng của chất khô hòa tan, Ms, được tính theo công thức của Schwartzberg (1976)
[94a]:
x s R g To2
(1.19)
M =
s

(Wwo − Wwb )L o t f

Trong đó thành phần nước liên kết trong thực phẩm, Wwb , được tính như sau:
(1.20)
Wwb = 0, 4Wp
Với, Wp là thành phần protein trong thực phẩm.
Thế phương trình (1.19) vào thành phần phân tử lượng của chất khô hòa tan trong phương
trình (1.18) ta thu được phương pháp đơn giản nhất dự đoán thành phần băng trong thực
phẩm (Miles 1974)[74]
 t 
(1.21)
W = (Wo − Wb ) 1 − f
ice

w

w



t 


trình biến đổi cân bằng pha có tỏa nhiệt. Trên cơ sở đó tác giả đã đề xuất công thức cải tiến
công thức (1.22) có dạng như sau:

13





Wice
1 
α
(1.23)
ω = o = 1 − o 1 −

ln(T
+
∆
−
T
)
Ww
Ww  1 + (α − 1)
f
ec 

ln(Tf + ∆ − T) 
Trong đó: α, ∆ là các hệ số thực nghiệm phụ thuộc từng loại thực phẩm, Tec- nhiệt độ
eutecti của thực phẩm.
Latyshev (1992)[65] giả thuyết rằng trong tất cả các trường hợp điểm eutecti của


Bien thien nhiet dung rieng theo nhiet do, t=[-40:40] C

Nhiet dung rieng, [kJ/kg.K]

100
Duoi diem dong bang
Tren diem dong bang
80

60

40

20

0
-40

-30

-20

-10

0

10

20

Khi không biết thành phần chi tiết cấu thành nên thực phẩm, Chen (1985)[33] đưa ra công
thức tính NDR của thực phẩm chưa kết đông một cách tương đối có dạng:

C u = 4190 − 2300.Ws − 628.Ws2

(1.27)
Các công thức (1.24÷1.27) khá đơn giản, dễ dàng để tính toán. Tuy nhiên lại không
phản ánh được sự phụ thuộc của NDR vào nhiệt độ. Do đó các công thức trên chỉ dùng để
tính toán định hướng trong các bài tính nhiệt, chứ không thể dùng để mô phỏng quá trình
cấp đông của thực phẩm nói chung và cá da trơn nói riêng.
Gupta (1990) [53] đề xuất phương trình xác định nhiệt dung riêngthực phẩm như là hàm
của nhiệt độ và thành phần nước trong dải nhiệt độ 303 ÷ 336K và thành phần nước 0,1%
÷80% như sau:
C = 2477 + 2356Ww – 3,79T
(1.28)
Mô hình toàn diện nhất xác định nhiệt dung riêng của các loại thực phẩm có kể đến
ảnh hưởng của nhiệt độ và các tất cả các thành phần trong thực phẩm được Choi và Okos
đưa ra năm 1986:[34]

C = Σ C i Wi

(1.29)
Trong đó Wi, Ci, là nhiệt dung riêng và thành phần khối lượng của thành phần thứ i
Phụ lục 1 trình bày nhiệt dung riêng của các thành phần trong thực phẩm là hàm phụ thuộc
nhiệt độ. Sự sai lệch giữa số liệu thí nghiệm và mô hình của Choi & Okos là do: (1) nhiệt
dung riêng của nước liên kết có sự sai khác rất lớn so với phần lớn lượng nước trong thực
phẩm, và (2) nhiệt dung riêng dôi ra do sự tác động qua lại của các pha thành phần.
Rahman (1993) [91] đã tính đến nhiệt dung riêng dôi ra để xác định C theo phương trình
(1.23), tuy nhiên việc này rất là khó thực hiện nếu chỉ dùng những phương pháp đơn giản.
Phương trình của Choi & Okos (1986)[34] chỉ ra rằng nhiệt dung riêng tăng khi nhiệt độ

(thường là -40oC)
Bằng mô hình đơn giản tương tự như của Schwartzberg (1976)[94a], Chen [33] đưa
ra mô hình xác định nhiệt dung riêng của thực phẩm kết đông bằng cách mở rộng phương
trình Siebel (1892)[95]:
Cf = 1550 + 1260.Ws +

Ws R g To2

(1.32)

Ms t 2

Nếu phân tử lượng của chất khô hòa tan chưa biết, ta có thể dùng công thức (1.19) để ước
lượng phân tử lượng của chất khô hòa tan, thế vào (1.32) ta được:
C f = 1550 + 1260.Ws −

(Wwo − Wwb )L o t f
t2

(1.33)

Ba mô hình dự đoán nhiệt dung riêng hiệu dụng (Chen 1985a[33],
Schwartzberg(1976) [94a], (1981)[94b]) đều cho kết quả tương tự và có sai số trung bình
tuyệt đối khoảng 20% và độ lệch tương đối lớn. Trong ba mô hình trên, mô hình của
Schwartzberg (1976)[94a] có sai số trung bình tuyệt đối thấp hơn so với hai phương trình
còn lại. Việc thực hiện tính toán theo mô hình của Schwartzberg (1981)[94b] khó khăn hơn
bởi vì nó dựa trên các giá trị nhiệt dung riêng của toàn bộ thực phẩm đông lạnh. Mô hình
dự đoán nhiệt dung riêng hiệu dụng của Chen (1985a)[33] là dễ sử dụng nhất, mặc dù có
sai số lớn nhất.
Mặt khác các công thức từ (1.30÷1.33) đều dựa trên công thức xác định hàm lượng


16


Lúc này, nước trong thực phẩm bắt đầu đóng băng và lượng băng hình thành bên
trong thực phẩm là hàm của nhiệt độ. Đồng thời trong biểu thức nhiệt dung riêng hiệu dụng
có xuất hiện thành phần của nhiệt ẩn chuyển pha của nước:
C e (T) = C(T) − Wwo .L(T)

dω
dT

Ce (T) = Ws .Cs (T) + Wwo (1 − ω).Cw (T) + Wwoω.Cice (T) − Wwo .L(T)

dω
dT

(1.36)

Trong đó:
ω - thành phần nước đóng băng được tính theo công thức (1.23)
L(T) - nhiệt ẩn đóng băng phụ thuộc nhiệt độ được tính như sau:
273,15

L(T) = L o −

∫
T

273,15

trình đo NDR dựa vào sự thay đổi nhiệt độ của mẫu đo và mẫu chuẩn kết hợp với phương
trình cân bằng nhiệt sẽ tính được nhiệt dung riêng của mẫu cần đo. Phương pháp xác định
NDR được thực hiện theo nhiều cách khác nhau gồm: Phương pháp hỗn hợp, phương
pháp nhiệt lượng kế so sánh, phương pháp buồng đoạn nhiệt, phương pháp
Differential scanning calorimeter (DSC)
Theo sơ đồ nguyên lý các phương pháp đo được trình bày trong phụ lục 2, các phương
pháp đo có ưu nhược điểm như sau:
Ưu điểm:
Các bộ thí nghiệm thiết kế đơn giản, ngoại trừ phương pháp DSC
Các thí nghiệm thực hiện nhanh, kích thước mẫu nhỏ
Nhược điểm:
Các thí nghiệm đo NDR cần phải có mẫu chuẩn
Tổn thất nhiệt lớn, tổn thất nhiệt trong quá trình đo không khống chế được
Không thể xác định được nhiệt dung riêng trong vùng kết đông của thực phẩm
Thiết bị đo cần phải có độ chính xác cao.
Phương pháp DSC có ưu điểm sau:
Các thí nghiệm thực hiện nhanh, số liệu chính xác.
Kích thước mẫu nhỏ
17


Số lượng mẫu cho thí nghiệm ít
Xác định được NDR của thực phẩm trong vùng kết đông
Nhược điểm:
Trong quá trình thí nghiệm cần phải có mẫu chuẩn để so sánh
Nhiệt độ mẫu đo cần phải đồng nhất
Thiết bị thí nghiệm phải hoàn toàn kín,tránh mất nước khi thực hiện các phép đo ở
nhiệt độ cao.
Chi phí thiết bị cao…
Kết luận: Trong các phương pháp thí nghiệm trình bày ở trên để xác định NDR


n

H = ∑ H W = ∑ ∫ C W dT
i =1

i

i

i =1

i

i

(1.39)

Từ mô hình nhiệt dung riêng của Chen (1985), enthalpy của thực phẩm chưa kết
đông được xác định bằng cách tích phân phương trình 1.21 theo nhiệt độ.
(1.40)
H = H f + (t − t f )(4190 − 2300 Ws − 628 Ws3 )
Trong đó Hf là enthalpy của thực phẩm ở điểm đóng băng, J/kg
Đối với thực phẩm dưới điểm kết đông, biểu thức toán học của enthalpy thu được
nhờ tích phân các biểu thức mô phỏng nhiệt dung riêng tương ứng. Tích phân phương trình
nhiệt dung riêng của Schwartzberg (1976)[94a] từ nhiệt độ tham chiếu Tr đến nhiệt độ của
thực phẩm T ta được phương trình: Schwartzberg

 R g T02


Ws R g T02

H = (t − t f )1550 + 1260.Ws −
M s .t.t f








(Wwo − Wbb )L 0 t f 
H = (t − t f )1550 + 1260 .Ws −

t.t f


Pham (1994) viết lại mô hình xác định enthalpy của Schwartzberg như sau:

H = A + CfT +

B
T

(1.42)
(1.43)

(1.44)



1.3.4

Hệ số dẫn nhiệt

1.3.4.1 Vai trò của hệ số dẫn nhiệt trong quá trình cấp đông
Hệ số dẫn nhiệt của thực phẩm ảnh hưởng bởi ba yếu tố: thành phần, cấu trúc và
điều kiện của quá trình. Thành phần nước cũng đóng vai trò rất quan trọng vì nước chiếm
một tỷ lệ rất lớn trong thực phẩm. Các yếu tố cấu trúc bao gồm độ rỗng, kích thước, hình
dáng và sự sắp xếp hay phân bố pha trong thực phẩm (khí, nước, băng và chất khô). Trong
đó độ rỗng đóng vai trò rất quan trọng vì hệ số dẫn nhiệt của chất khí nhỏ hơn của chất
lỏng và rắn tới hàng chục thậm chí cả trăm lần. Các yếu tố quá trình bao gồm nhiệt độ, áp
suất và điều kiện truyền nhiệt. Khác với định nghĩa của hệ số dẫn nhiệt theo các lý thuyết
truyền nhiệt kinh điển, hệ số dẫn nhiệt của thực phẩm có thể thay đổi theo phương vuông
góc hay song song với các thớ trong thực phẩm, do đó việc dự đoán và đo hệ số dẫn nhiệt
của thực phẩm là rất phức tạp và hệ số đo được ở đây là hệ số dẫn nhiệt hiệu dụng.

1.3.4.2 Mô hình toán dự đoán hệ số dẫn nhiệt
Maxwell (1904) [68] là người đi tiên phong trong lĩnh vực nghiên cứu hệ số dẫn
nhiệt của hỗn hợp hai pha. Phương trình của Maxwell dựa trên lý thuyết điện thế được viết
như sau:
19


 λ + 2λc − 2ε(λc − λd ) 
(1.46)
λ = λc  d

λ
+

VC + Vd
Kopelman (1966) [64] đã nghiên cứu về mô hình truyền nhiệt trong các loại thực
phẩm có thớ và tìm ra sự khác biệt khi nhiệt truyền dọc thớ và ngang thớ thịt. Đối với hệ
đẳng hướng, hệ hai thành phần được cấu thành từ pha liên tục và pha gián đoạn, hệ số dẫn
nhiệt độc lập với hướng dòng nhiệt, Kopelman đưa ra biểu thức hệ số dẫn nhiệt sau:
a=



1 − ε2/3
λ = λc 
2/3
1/3 
1 − ε (1 − ε ) 

Trong đó :
ε là thành phần thể tích của pha gián đoạn.
Theo (1.47), hệ số dẫn nhiệt của pha liên tục được cho là lớn hơn rất nhiều so với
hệ số dẫn nhiệt của pha gián đoạn. Tuy nhiên, nếu hệ số dẫn nhiệt của pha gián đoạn lớn
hơn nhiều hệ số dẫn nhiệt của pha liên tục thì công thức sau được sử dụng:


1− ξ
λ = λc 
1/3 
1 − ξ(1 − ε ) 

(1.48)

λd


(1.50)

ϕ = ε1/2 (1 − λd / λc )
Công thức của Kopeman đưa ra có kể đến hướng truyền của dòng nhiệt. Ứng với
hai phương của dòng nhiệt, ta có thể giải được bài toán dẫn nhiệt một chiều. Trong thực tế
20


thì dòng nhiệt là vô hướng và cần phải sử dụng đến hệ số dẫn nhiệt hiệu dụng thì từ hai mô
hình theo hai phương của dòng nhiệt lại chưa đưa ra được.
Levy (1981) [67] đưa ra dạng khác kết hợp cả công thức Maxwell và Eucken cho
hệ hai thành phần như sau:
λ [ (2 + Λ ) + 2( Λ − 1)F1 ]
(1.51)
λ= 2
(2 + Λ ) − ( Λ − 1)F1

Trong đó:

Λ = λ1 / λ 2 - tỷ số hệ số truyền nhiệt của thành phần 1 và 2
0,5
2
 2


  2
 8R 1  
F1 = 0,5  − 1 + 2R 1  −  − 1 + 2R1  −
 

2



Khi thực phẩm bao gồm nhiều hơn hai pha riêng biệt, phương pháp đề ra trên đây
để dự đoán hệ số dẫn nhiệt được ứng dụng để tính hệ số dẫn nhiệt của thực phẩm. Ví dụ,
trong trường hợp thực phẩm kết đông, hệ số dẫn nhiệt của băng và nước lỏng được tính
toán trước bằng phương pháp trên. Hệ số dẫn nhiệt thu được của hỗn hợp băng nước này
lần lượt được kết hợp với hệ số dẫn nhiệt của từng thành phần trong thực phẩm để xác định
được hệ số dẫn nhiệt của toàn khối thực phẩm. Phương pháp nêu trên hiện nay được đánh
giá là cho sai số khi tính toán là thấp nhất.
Một số các nhà nghiên cứu đưa ra mô hình dẫn nhiệt song song và dọc theo thớ
thực phẩm dựa trên mô hình tương tự điện trở (Murakami và Okos 1989) [76]. Mô hình
song song, hệ số dẫn nhiệt tổng là tổng của tích các hệ số dẫn nhiệt thành phần với thành
phần thể tích:
(1.55)
λ pa = ΣVi λ i
Trong đó Vi là thành phần thể tích của thành phần thứ i,
Wi / ρi
Vi =
Σ(Wi / ρi )
Ở mô hình nối tiếp, dòng nhiệt vuông góc với thớ, hệ số dẫn nhiệt bằng nghịch đảo
của tổng các thương giữa thành phần thể tích và hệ số dẫn nhiệt thành phần:
1
(1.56)
λ =
se

Σ(Vi / λ i )



λ = −q

∆x
∆T

Từ đ ó

(1.57)

(1.58)
Bằng cách tạo ra dòng nhiệt ổn định một chiều qua vật thể có điều kiện biên loại 1
ổn định, hệ số dẫn nhiệt được xác định khi đo được mật độ dòng nhiệt dẫn qua vật và nhiệt
độ trên hai mặt ngoài vật. Dòng nhiệt dẫn qua vật thí nghiệm được xem là công suất nhiệt
của dòng điện đốt nóng. Điều quan trọng là chọn dạng vật đo và thiết bị đo sao cho phù
hợp để dòng nhiệt là một chiều.
Phương pháp ổn định có ưu điểm: xác định được hệ số dẫn nhiệt có độ chính xác
tương đối cao trong phạm vi nhiệt độ lớn, quá trình tính toán và xử lý đơn giản. Tuy nhiên
phương pháp này có một số nhược điểm:
Thời gian thí nghiệm rất lâu: vài chục giờ đến vài ngày.
Hệ thống đo cồng kềnh,mẫu cần đo bắt buộc phải gia công.
Chỉ xác định được λ trung bình mà không xác định được λ tức thời.
Vật liệu cần đo phải khô, vì khó có thể tạo được chế độ ổn định với vật ẩm.
Không thể xác định λ của vật có trị số dẫn nhiệt rất nhỏ (λ
mặc dù được pháp triển dưới rất nhiều dạng khác nhau. Phổ biến rộng rãi hơn là phương
pháp dây đốt do Schleidermacher đề xuất 1880. Đây là phương pháp tiền thân của tất cả
các phương pháp được biết đến dưới tên gọi là phương pháp”que thăm”sau này. Như đã đề
cập ở trên phương pháp này tuy có thể dùng để đo hệ số dẫn nhiệt của thực phẩm cấp đông,
nhưng quá trình thực nghiệm của phương pháp này khá phức tạp, rất khó xác định chính
xác lượng nhiệt truyền qua các lớp vật liệu do không xác định được thật chính xác nhiệt
lượng tích tại nguồn.
Vì những lý do trên trong khuôn khổ nội dung nghiên cứu của luận án này tác giả đề
xuất áp dụng phương pháp xung nhiệt để xác định bằng thực nghiệm hệ số dẫn nhiệt hiệu
dụng λ(T) và hệ số dẫn nhiệt độ a(T).

1.3.5

Khối lượng riêng

Khối lượng riêng của thực phẩm phụ thuộc chủ yếu vào thành phần thực phẩm và
nhiệt độ. Thực phẩm được coi là một hệ nhiều pha, trong đó khối lượng và thể tích là
không đổi trong các quá trình hỗn hợp. Choi & Okos (1986)[34] để xuất phương trình dưới
đây xác định khối lượng riêng của thực phẩm khi đã biết thành phần khối lượng của chúng:
1
(1.59)
ρ=
Wi
∑i ρ
i
Tuy nhiên, phương trình trên không kể đến pha khí cũng như sự tương tác giữa các
pha, do đó ứng dụng của nó bị hạn chế. Trong trường hợp sản phẩm rỗng, với độ rỗng ε,
khối lượng riêng được tính theo công thức:
(1 − ε )
(1.60)

trường nhiệt độ sớm hơn so với vật liệu có hệ số dẫn nhiệt độ bé, do đó ta có mối quan hệ
sau:

a
∆t
=
a
∆t
1

2

2

1

(1.62)

Trong đó: a1,a2: hệ số dẫn nhiệt độ của vật liệu 1 và vật liệu 2[m2/s]
∆t1,∆t2: Thời gian cần thiết để đạt được cùng một trường nhiệt độ[S]
1.3.6.2 Phương pháp xác định
Về phương pháp xác định hệ số dẫn nhiệt độ được chia ra thành 2 phương pháp:
phương pháp xác định gián tiếp và phương pháp xác định trực tiếp
a.
Phương pháp xác định gián tiếp
Theo phương pháp này, hệ số dẫn nhiệt độ được xác định thông qua 3 thông số
nhiệt vật lý khác là: hệ số dẫn nhiệt λ, nhiệt dung riêng C và khối lượng riêng, được tính
như sau:

a=


(1.64)

Trong đó β: tham số không thứ nguyên;
r
β=
2 (a × τ)

(1.65)

Với:
r: khoảng cách giữa hai que thăm [m]
a: hệ số dẫn nhiệt độ [m2/s]
t: Thời gian [s]
Theo Nix và cộng sự (1967), lời giải của phương trình (1.64) như sau:
q
β2
β4
 Ce

∆T =
−
− ln(β) +
−
+ ...
(1.66)

2× π× λ  2
2 × 1! 4 × 2!
