Kỹ thuật tính toán trong công nghệ hóa học PDF

Kỹ thuật tính toán ừong công nghệ hóa học
Mục Lục
Lời nói đầu............................................................................................................................... 3
Phần I: Cơ sở về Matlab..........................................................................................................4
Chương I:Bắt đầu VỔ4 Matlab................................................................................................ 5
1.1 Nhập dữ liệu qua dòng lệnh............................................................................................. 5
1.2 Sử dụng help trực tuyến....................................................................................................6
1.3 Đường dẫn.........................................................................................................................7
1.4 Lưu và tải dữ liệu............................................................................................................. 7
Chương 2: Các cấu trúc cơ bản và biến................................................................................ 8
2.1 tính toán với Matlab......................................................................................................... 8
2.2 Giới thiệu dạng dữ liệu.....................................................................................................9
Chương 3: Tính toán vói ma trận và vec tơ .......................................................................10
3.1 Vec tơ............................... ............................................................................................... 10
3.1.2 Vec tơ cột và sự chuyển vị................................................................................. 12
3.1.3 Nhăn, chia và mũ của vec tơ .............................................................................12
3.2 M aừ ận............................................................................................................................ 14
3.2.1 Các ma frận đặc biệt..........................................................................................15
3.2.2 Xây dựng ma ừận và cách trích các phần tử của ma trận............................. 16
3.2.3 Các phép toán với ma trận................................................................................19
Chương 4: Đồ họa...................................................................................................................22
4.1 Đồ thị đơn giản............................................................................................................... 22
4.2 Vẽ đường........................................................................................................................ 23
4.3 Vẽ mặt............................................................................................................................. 26
Chương 5: Các luồng điều khiển..........................................................................................30
5.1 Các toán tử logic............................................................................................................ 30
5.2 Lệnh find......................................................................................................................... 30
5.3 Cậu lệnh if.......................................................................................................................32
5.4 Cấu trúc lặp.....................................................................................................................34
Chương 6: Phương pháp sổ...................................................................................................35
6.1 Đường cong hồi quy....................................................................................................... 35
6.2 Phép nội suy....................................................................................................................35
6.3 Giá trị của hàm số........................................................................................................... 36
6.3.1 Hàm 'inline'....................................................................................................... 37
6.4 Phép tính tích phân vạ tích phân....................................................................................38
6.5 Tính toán số học và cấu trúc điều khiển........................................................................38
Chương 7: Viết chương trình trong Matlab.......................................................................39
7.1 Kịch bản và hàm m-fĩle................................................................................................. 39
7.2 Hàm m-file......................................................................................................................39
7.2.1 Những biến đặc biệt của hàm s ố ......................................................................40
7.2.2 Biến địa phương và biến toàn cục.................................................................... 41
7.2.3 Tính giá trị hàm số cách gián tiếp.................................................................... 41
1
Chương 8: Văn bân.................................................................................................................43

8.1 Chuỗi kí tự ...................................................................................................................... 43
8.2 Vào - ra dữ liệu .............................................................................................................. 44
8.2.1 Vào ra dữ ỉiệu từ bàn p h ím ..............................................................................44
8.2.2 Đưa dữ ỉiệu ra màn hình.................................................................................. 44
8.2.3 Xứ lífile văn bản...............................................................................................46
Phần II: ứng dụng của Matlab............................................................................................ 51
trong công nghệ hóa học........................................................................................................ 51
Chương 1: Dẩn nhiệt và đối lưu........................................................................................... 52
1.1. Dẩn nhiệt................ ....................................................................................................... 52
1.1.1. Tóm tắt lí thuyết.............................................................................................. 52
1.1.2. Vỉ dụ................................................................................................................. 53
1.2. Đối lưu.......................................................................................................................... 56
1.2.1. Tóm tắt lí thuyết............................................................................................... 56
ỉ.2.2. Vỉ dụ.................................................................................................................. 59
Chương 2: Bức xạ nhiệt và truyền nhiệt............................................................................. 64
2.1. Trao đổi nhiệt bức xạ giữa các vật............................................................................... 64
2.1.1. Tóm tắt lí thuyết............................................................................................... 64
2.1.2. V ỉđụ.................................................................................................................64
2.2. Truyền nhiệt và thiết bị trao đổi nhiệt..........................................................................67
2.2.1 Tóm tắt ỉí thuyết................................................................................................ 67
2.2.2. Vỉ d ụ ................................................................................................................. 68
Chương 3: Kỹ thuật tách chất.............................................................................................. 74
3.1 Cân bằng lỏng h ơ i.......................................................................................................... 74
3.1.1 Tóm tắt lỷ thuyết............................................................................................... 74
3.1.2 Ví dụ.... ..............................................................................................................75
3.1.3. Bài tập.............................................................................................................. 79
3.2 Phương pháp McCabe Thiele........................................................................................ 80
3.2.1 Tóm tắt lý thuyết............................................................................................... 80
3.2.2 Vỉ dụ.... ..............................................................................................................83
Chương 4: Kỹ thuật phản ứng............................................................................................. 89
4.1 Nhiệt động học................................................................................................................89
4.1.1 Tóm tắt ỉỷ thuyết............................................................................................... 89
4.1.2 Vi dụ.... ..............................................................................................................91
4.2 Các kỹ thuật tiến hành phản ứng....................................................................................94
4.2.1 Tóm tắt lý thuyết............................................................................................... 94
4.2.2 Vỉ dụ .... ..............................................................................................................95
4.2.3BÒÌ tập................................................................................................................ 99
Phụ lục.................................................................................................................................... 100
Tài liệu tham khảo................................................................................................................125
2
Lời nói đầu
Các bài toán công nghệ hóa học trong thực tế là rất phức tạp, chúng ta cần có
công cụ hỗ ừợ để giải nó. Thực tế có rất nhiều phần mền, tuy nhiên, ừong cuốn sách
này, chúng tôi muốn giới thiệu đến các sinh viên nghành công nghệ hóa một công cụ
rất mạnh, đó là matlab .
Mat lab là một công cụ toán học rất mạnh, ngoài ra nó còn hỗ trợ ngôn ngữ lập
trình bậc 4 vói cấu trúc đơn giản, gần gũi, dễ tiếp cận, có thể lỉnh hoạt giải quyết các
bài toán thực tế để đạt được yêu cầu mong muốn. Do đó, việc sử dụng nó làm công cụ
để giải quyết các bài toán kĩ thuật là rất tố t.
Nội dung cuốn sách bao gồm:
• Phần I: Cơ sở về Matlab nhằm giới thiệu cho bạn đọc sơ lược về Matlab
• Phần II: ứ n g dụng Matlab ừong công nghệ hóa học gồm có các chương:
Dẩn nhiệt và đối lưu; bức xạ và truyền nhiệt; kỹ thuật tách chất; kỹ thuật
phản ứ n g .
Trong từng chương có tóm tắt cơ sở lí thuyết, kèm theo mỗi chương đó thì có các
ví dụ và bài tập liên quan đi kèm, mỗi một ví dụ minh họa được giải bằng tay sau đó
được giải lại bằng matlab một cách sinh động và dễ hiểu, có tính tổng quát cao .
Đưa matlab vào công nghệ hóa là mục đích của chúng tôi. Tuy nhiên, do được
biên soạn lần đầu nên không thể ừánh được những sai sót, rất mong nhận được đóng
góp của bạn đọc để cuốn sách trở lên gần gũi hơn với các sinh viên công nghệ hóa
cũng như những người quan tâm đến lmh vực n à y .
3
Phần I: Cơ sở về Matlab
Matlab là công cụ toán học, trước hết là các phép tính khoa học, ngoài ra nó còn
cho phép vẽ các đồ thị, tính toán ma trận, làm việc với các đa thức và các hàm tích
phân. Đồng thời,chúng ta có thể lập trình các chương trình tính toán một cách linh hoạt
theo nhu cầu của mình. Đặc biệt Matlab tỏ ra rất manh mẽ trong việc thiết kế và tính
toán vector và ma trận.
Khi khởi chạy Matlab thì xuất hiện một hoặc nhiều cửa sổ trên màn hình, trong đó
quan trọng nhất là Command window, nó cho phép nhập lệnh và hiển thị kết quả.
4
Chương I:Bắt đầu với Matỉab
1.1 Nhập dữ liệu qua dòng lệnh

*) Các lệnh sẽ được thực hiện ngay sau khi ấn phún Enter. Ket quả sẽ được hiển
thị trên màn hình nếu muốn. Tuy nhiên, nó chỉ thực hiện khi các dòng lệnh là hợp lệ.
Có thể tìm hiểu thêm từ menu Help.
Ví d ụ :
» 3 + 7.5
» 18/4
» 3 * 7
Chú ý : các dấu cách không quan ừọng trong Matlab
*) Kết quả tính toán sê được tự động lưu vào biến ans
Ví d ụ :
» 14/4
ans =
3.5000
» ansA(-6)
ans -
5.4399e-04
Chúng ta có thể định nghĩa biến cho riêng mình. Ví dụ:
» a = 14/4
a
3.5000
» b = aA(-6)
b=
5.4399e-04
Dữ liệu sẽ được lưu vào biến a và b.
*) Khi kết thúc dòng lệnh kết thúc bằng dấu thì sẽ không hiển thị kết quả ra
màn hỉnh, có thể kiểm chứng bằng hai biểu thức sau :
» 3 + 7.5
» 3 + 7.5;
*) Có thể thực thi nhiều lệnh cùng một lúc, các câu lệnh ngăn cách nhau bởi dấu
phẩy hoặc dấu chấm phẩy.
>> sìnịpì/4'), cos(pi); sin(O)
5
ans
0.707ì
ans
0
Chú ý : cos(pi) không được in ra màn hỉnh
*) Mặc định Matlab chỉ hiển thị 5 chữ số. Do đó, muốn có kết quả chính xác hơn
thì ta dùng lệnh form at long thì có thể hiển thị được tới 15 chữ số
» 312/56
ans =
5.5714
» format long
» 312/56
ans
5.57142857142857
Có thể tìm hiểu thêm:
>> help format
*) Khi dòng quá dài có thể dùng dấu để nối câu xuống dòng sau
>> sin(l) + sin(2) - sin(3) + sin(4) - sin(5) + sin(6) - ...
sin(8) + sin(9) - sin(10) + sin (ll) - sin(12)
ans
1.0357
*) Dòng kí tự đặt sau dấu “%” chỉ có tác dụng làm dòng chú thích
>> sin(pi) % dùng để tính giá trị của sin (pi)
1.2 Sử dụng help trực tuyến

Chúng ta có thể dùng cấu trúc
>> help <cầu ỉệnh>
Vd: help ops
Hoặc
>> lookfor <câu ỉệnh>
Vd : lookfor inverse
6
1.3 Đường dẫn

Trong Matlab, các lệnh hay các chương trình được chứa trong m-file, cái mà được
biên soạn từ file văn bản nhưng khi lưu vói đuôi mở rộng là “.m”. M-file chi chạy
được kbi nó nằm ở thư mục current working directory
1.4 Lưu và tải dữ liệu

Rất dễ dàng lưu hoặc gọi 1 biến nào đó, ta chỉ can click vào File- menu sau đó
chọn Save Workspace as... hoặc Load Workspace...
V d : khi save
>> s i = sin(pi/4);
» c l = cos(pi/4); c2 = COS(pi/2);
» str= ’hello world’; % đây là 1 chuỗi
» save % lưu toàn bộ biến với định dạng matlab.mat
» save data % lĩtu các biến với đpth dạng nhị phân data.mat
» save numdata si, cl % lưu các biến số sl và c l tới numdata.mat
» save strdata sfr % hru biển chuỗi Sừ vào strdata.mat
» save allcos.dat c* -ascii % saves cl,c2 in 8-digit ascii format to
allcos.dat
Còn khi load
» load % tải tẩt cả các biến từfile matlab.mat
» load data s i cl % chi tải riêng từng biển trong file data.mat
» load Sừ- data % tải tất cả các biến từ file strdata.mat
7
Chương 2: Các cấu trúc cơ bản và biến
2.1 Tính toán với Matlab

Có 3 loại biến số học: số nguyên, số thực, số phức. Matlab còn hỗ trợ thêm
những loại không phải là số: Inf (giá trị vô cực, ví dụ : 1/0 sẽ cho kết quả như vậy)
hoặc NaN (không phải là số, là kết quả của các phép tính như 0/0 hoặc vô cùng - vô
cùng)
Chúng ta có thể viết các biểu thức trên dòng lệnh của Matlab :
» (23*17)/7
Kết quả sẽ là :
ans
55.857!
Matlab có 6 toán tử cơ bản là \ , /(chia phải hoặc ừái) và A(hàm mũ)
Chú ý sự khác nhau của toán tử chia trái và chia phải :
>> 19/3 % thực hiện phép chia: 19/3
ans
6.3333
» 19\3, 3/19 % thực hiện phép chia: 3/19
ans =
0.1579
ans
0.1579
Những hàm cơ bản về lượng giác, mũ ... có thể tìm hiểu
>> help elfim
Bài tập:
Xác định giá trị các biểu thức dưới đây bằng tay và bằng Matlab để so sánh. Chú
ý tới sự khác nhau giữa chia trái và phải, và dùng lệnh help để tìm hiểu về các hàm:
round, floor, ce il...
•2/2-3
• 8 - 5\4
• 8 - (5\4)
• 7 - 5 -4\9
• 6 - 2/5 + 7A2 - 1
8
• 10/2\5 - 3 + 2 - 4
• 3A2/4
• 3A2A3
• 2 + round (6/9 + 3 - 2)12
• 2 + floor (6/9 + 3 - 2)/2
• 2 + ceil (6/9 + 3 - 2)12
• X = pi/3, x = x - l , x = x + 5, x = abs(x)/x
2.2 Giới thiệu dạng dữ liệu

Mặc dù tất cả các tính toán số trong Matlab đều được tính toán với độ chính xác
kép (double precision), nhưng khuôn dạng của dữ liệu đưa ra có thể địnhdạng lại nhờ
các lệnh định dạng của Matlab. Các biến ngầm định cũng nhưcác biến của người sử
dụng định nghĩa đều có thể đưa ra nhiều định dạng khác nhau. Đinh dạng được chọn
nhờ sử dụng lệnh form at tiếp đến là chỉ định long hay short và cuối cùng có thể là
dạng biểu diễn của dữ liệu, cụ thể là e đối với dạng hàm mũ, f đối dạng dấu phẩy
động(ngầm định là f). Một khuôn dạng ít dùng hơn là rat để đưa dữ liệu ra dưới dạng
phân số (chính xác hoặc gần đúng ).
Ví dụ:
>> 4*atan(l)
ans =
3.1416
» format longe
» ans
ans =
3.141592653589793e+000
» format lo n g f
» ans
ans
3.14159265358979
» format short e
» ans
ans
3.1416c ■000
» format rat
» ans
ans
355/113
9
Chương 3: Tính toán vói matrận và vec tơ
Nên tảng của Matlab chính là ma trận (hay mảng).Trường họp đặc biệt:
1 ma trận 1X1: là 1 số vô hướng;
1 ma trận chỉ có 1 dòng hoặc 1 cột: là 1 vec tơ
3.1 Vec tơ
Vec tơ hàng là 1 danh sách các số được ngăn cách nhau bởi dấu cách hoặc dấu
phẩy. Số lượng phần tử của 1 vec tơ được gọi là length của vec tơ. Các phần tử nhập
vào phải đặt ừong dấu “[ ]”. Vd:
>> V = [-1 sm(3) 7]
V=
-1.0000 0.1411 7.0000
» length(v)
ans
3
Một số phép toán có thể thực hiện với vec tơ như nhân vô hướng, cộng hoặc trừ
từ các vec tơ khác có số phần tử tương tự, hoặc giữa 1 số với vec tơ. Mọi phép toán
đều được thực hiện trcn từng phần tử. Các vec tơ cũng có thể được xây dựng từ các
phần tử hữu hạn.
» V = [-1 2 7J; w = [2 3 4];
» z - V+ w % cộng từng phần tử
1 511
» vv = V + 2 % cộng 2 vào tất cả các phần tử vào vec tơ
vv =
149
» t = [2*v, -wj
ans -
-2 4 14 -2 -3 -4
Đồng thời 1 phần tử có thể được trích dẫn ra, hoặc chinh sửa:
» v(2) = -1 % sửa phần tử thứ 2 của V
-1 -1 7
» w(2) % hiển thị phần tử thứ 2 của w
10
ans -
3
3.1.1 Kỉ tự hai chẩm và ừích các phần tử của vec tơ

Dấu hai chấm là phím tắt quan ừọng dùng để tạo ra vec tơ hàng:
» 2:5
ans =
2 3 4 5
» -2:3
ans
2 ì 0 1 2 3
Dạng tổng quát “first:step:last”, khởi tạo vec tơ dòng bắt đầu từ phần tử first, kết
thúc tại phần tử last và với độ dài bước là step. Nếu không chỉ ra step, thì giá trị mặc
định nó là bằng một.
» 0.2:0.5:2.4
ans
0.2000 0.7000 1.2000 1.7000 2.2000
» -3:3:10
ans =
-3 0 3 6 9
» 1.5:-0.5:-0.5 % cũng có thể dùng bước nhảy âm
% -1:2:6 => -1 1 3 5
» r(3:6) % lẩy ra các phần từ của r cổ vị tri từ 3 đến 6

ans
3 5 2 3
» r(l:2:5) %lấy ra các phần tử của r có vị trí 1,3 và 5
ans =
- 1 3 2
11
>> r(5:-l:2) % chuyện gì sẽ xảy ra với trường hợp này?
3.1.2 Vec tơ cột và sự chuyầĩ vị

Để tạo 1 vec tơ cột, chúng ta phải tách các phần tử bởi phím ‘enter’ hoặc dấu
chấm phẩy:
» z = [1
7
7];
z=
1
7
7
» u = [-1; 3; 5]
u=
-1
3
5
Các phép toán cũng được thực hiện 1 cách tuơng tự, và cũng có thể chuyển dạng
thành vec tơ hàng:
» u ’ % u là vec tơ hàng cỏn u ’ là vec tơ cột
ans =
- 1 3 5
3.1.3 Nhân, chia và mũ của vec tơ

Chúng ta có thể tính tích vô hướng của 2 vec tơ có cùng độ dài:
» u = [-1; 3; 5] % 1 vec tơ cột
» V = [-1; 2; 7] % 1 vec tơ cột
» u * v % chúng ta không thể nhân vec tơ cột với vecto cột
??? Error using ——> *
Inner maừix dimensions must agree.
» u’*V % đây mới là tích vô hướng
ans =
Một cách khác để tính tích vô hướng là sử dụng dấu

>> u .* V %nhân từng phần tử cho nhau
1
12
6
35
» sum(u. *v) % cách lẩy tích vô hướng
ans
42
» z = [4 3 1]; % z l à vec tơ hàng
» sum(u *z) % tích vô hướng
ans
Kì
Trong toán học, chứng ta không chia vec tơ cho vec tơ. Tuy nhiên, ừong Matlab
chúng ta có thể chia từng phần tử cho nhau của các vec tơ có cùng kích cỡ và cùng
loại:
>> Jt = 2:2:10
X
2 4 6 8 10
» y - 6:10
y =
67 8910
» x./y
ans
0.3333 0.5714 0.7500 0.8889 1.0000
» z = -1:3
z—
-1 0 1 2 3
» x./z % chia 4/0, kết quả ra In f
Warning: Divide by zero.
ans
-2.0000 In f 6.0000 4.0000 3.3333
» z./z % chia 0/0, kết quả ra NaN
13
ans
1 NaN 1 1 1
Toán tử có thể thực hiện bởi 1 đại lượng vô hướng với 1 vector:
» x=l:5; 2/x % không thể thực hiện
??? Error using ==> /

Matrix dimensions must agree.
» 2./x % nhimg thế này thì được!
ans
2.0000 1.0000 0.6667 0.5000 0.4000
3.2 M a trận
Một ma trận
'l 2 3Ì
A= 4 5
7 8
Được định nghĩa trong Matlab như sau:
»A = [l 2 3; 4 5 6; 7 8 9] % nhập theo hàng
A =
1 2 3
4 5 6
7 8 9
Hoặc cũng có thể:
» A 2 = [1:4; -1:2:5]
A2 =
1 2 3 4
-1 1 3 5
» A3 = [1 3
-4 7]
A3 =
1 3
-4 7
Phép toán chuyển vị ma ừận cũng có thể thực hiện 1 cách dễ dàng:
»A 2
A2 =
1 2 3 4
-1 1 3 5
» A 2’ % chuyển vị A2
ans =
14
/ -1
2 1
3 3
45
» size(A2) % lấy cỡ của A2
ans
2 4
» size(A2’)
ans =
4 2
3.2.1 Các ma trận đặc biệt

» E = [] % tạo lập 1 ma trận rỗng!
E =
[]
» size(E)
ans
0 0
» I = eye(3); % ma trận đơn vị 3x3
/ =
16
1 0 0
0 1 0
0 0 1
» X - [2; -1; 7J; I*x % I is such thatfor any 3-by-l X holds I*x = X
am -
2
-1
7
» r - [1 3 -2J; R = diag(r) % tạo lập 1 ma trận chéo với r nằm trên đường
chéo
R =
1 0 0
0 3 0
15
0 0 - 2
» A = [1 2 3; 4 5 6; 7 8 9];
» diag(A) % ừích dẫn đường chéo của ma trận A
ans
ì
5
9
» B = ones(3,2)
B =
1 1
1 1
1 1
» c = zeros (size(C‘) ) % một ma trận toàn phần tử toàn là sổ 0 với
kích thước bằng ma ừận C ’
c=
0 0 0
0 0 0
» D = rand(2,3) % tạo 1 ma trận với các phần tử là các so ngẫu nhiên
D=
0.0227 0.9101 0.9222
0.0299 0.0640 0.3309
» V = linspace(l,2,4) % mộtvec tơ được sử dụng như 1 ma ừận
V=
1.0000 1.3333 1.6667 2.0000
3.2.2 Xây dựng ma trận và cách trích các phần tử của ma trận
Người ta thường xây dựng một ma trận lớn hơn từ những cái nhỏ hơn:
» x = [ 4 ; - l ] , y = [-l 3]
X =
4
-1
y =
-1 3
»X= [xy7 % Xg ồm có c ộ t x v à y ’
16
x=
4-1
-1 3
» T - [-1 3 4; 4 5 6J; t = l:3;
» T - [T; t] % thêm hàng cho T, hàng đó là vec tơ t
T=
- 1 3 4
4 5 6
1 2 3
» G = [1 5; 4 5; 0 2];% G là m a trận 2*3
» T2 = [T GJ % ghép 2 ma trận
T2 =
-1 3 4 1 5
4 5 6 4 5
1 2 3 0 2
» T ỉ = [T; Gones(3,ỉ)J % G l à 3*2, Tlà 3*3
T3 =
- 1 3 4
4 5 6
1 2 3
1 5 1
4 5 1
0 2 1
» T 3 = [T; G ’] ;
» [G ‘ diag(5:6); ones(3,2) TJ % chúng ta cso thể nối nhiều ma trận
am
4 0 5 0
5 2 0 6
1 - 1 3 4
1 4 5 6
1 1 2 3
17
Một phần có thể được trích tù 1 ma trận tương tự như cách lấy của vectơ. Mỗi
phần tử trong ma trận được mang chỉ số bởi vị trí của nó trong hàng và cột. Ví dụ:
» A = [1:3; 4:6; 7:9]
A =
123
4 56
7 89
» Ả (1 ,2 ),A (2 ,3 ),A (3 ,1 )
ans -
2
ans -
6
arts
7
» A(4,3) % this is not possible: A is a 3-by-3 matrixl

??? Index exceeds matrix dimensions.
» A(2,3) - A(2,3) + 2*A(1,1) % sửa giá trị của phần tửA(2,3)
A =
1 2 3
4 5 8
7 8 9
Ta có thể dễ dàng mở rộng ma trận:
>> A(5,2) = 5 % gán 5 vào vị ừi (5,2);
Á =
1 2 3
4 5 8
7 8 9
0 0 0
0 5 0
» A(4,:) = [2,1, 2]; % gán vec tơ [2, 1, 2] cho hàng thứ 4 của A
» A(5,fl,3J) = [4, 4]; % gán : A(5,l) = 4 và A(5,3) = 4
» A
Lấy các phần tử khác của A:
>> A(3,:) % ừich ra hàng thứ 3 của A
18
ans -
7 8 9
» Ả(:,2) % trích ra cột thứ 2 của A
ans
2
5
8
1
5
» A(1:2,:) % trích ra hàng 1 và 2 của A
ans =
1 2 3
4 5 8
» A([2,5],l:2) % trích ra các phần của A
ans -
4 5
4 5
Khái niệm về ma trận rỗng [] rất hữu ích trong Matlab. 1 số hàng hoặc cột của ma
trận có thể được loại bỏ bằng cách gán nó bằng 1 ma trận rỗng:
» c = [l 2 3 4; 5 6 7 8;1 1 1 1];
» D = C; D(:,2) = [] % copy c vào D, xóa cột thứ 2 của D
% xóa hàng 1 và hàng 3 của c
3.2.3 Các phép toán với ma trận

Các phép toán thực hiện với ma trận cũng tuơng tự với vec tơ:
>> B = [1 -1 3; 4 0 7]
B =
1 -1 3
4 0 7
» B2 - [1 2; 5 1; 5 6];
» B = B + B2’ % cộng 2 ma trận. Tại sao phải thay B2 ’ cho B2?
B =
2 4 8
6 1 13
19
> > B-2 % ừ ừ các phần tử của B cho 2

ans -
0 2 6
4 -1 11
» ans - B./4 % chia các phần tử của B cho 4
ans
0.5000 1.0000 2.0000
1.5000 0.2500 3.2500
» 4/B % đây là điầi không thế
??? Error using ==> /
Matrix dimensions must agree.
» 4./B % nhưng đây lại có thể ; tương đương với: 4. *ones(size(B')) ./ B
ans -
2.0000 1.0000 0.5000
0.6667 4.0000 0.3077
» c - [1 -1 4; 7 0 -1];
» B .* c % nhân tùng phần tử
ans
2 -4 32
42 0 -13
» ans. ^3 - 2 % thực hiện cho tất cả các phần tử
ans -
6 -66 327 66
74086 -2 -2199 20
» a n s ./ B.A2
ans
0.7500 -1.0312 63.9961
342.9907 -2.0000 -1.0009
» r - [1 3 -2J; r * B2
ans
ù -7
Khi thực hiện nhân ma ừận:
» b = [13 -2];
» B = [1 -1 3; 4 0 7]
20
B =
1 -1 3
4 0 7
» b *B % không thể thực hiện: b là ma trận 1*3 và B là 2*3
??? Errorusmg ==> *
Inner matrìx dimensions must agree.
» b *B’
am
-8 -10
» B ’ *ones(2,l)
am =
» c = [3 1; 1 -3];
»C *B
ans
7 -3 16
-11 -1 -18
» C.A3% thực hiện ừên từng phần tử
ans
27 1
1 -27
» CA3 % tương đương với c*c*c
ans
30 10
10 -30
» ones(3,4)./4 * diag(l:4)
ans =
0.2500 0.5000 0.7500 1.0000
0.2500 0.5000 0.7500 1.0000
0.2500 0.5000 0.7500 1.0000
21
Kỹ thuật tính toán trong công nghệ hóa học
Chương 4: Đồ họa
4.1 Đồ thị đom giản

Đơn gỉẳn nhất, đồ thị có thể được vẽ nhờ nối các điểm được đánh dấu trên mặt
phẳng tọa độ Đề các. Ví dụ:
>> X = 0:10;
» y = 2.*x;% ta được y = [1 2 4 8 1 6 32 64 128 256 512 1024]
» plot(x,y) % vẽ
» semilogy(x,y)
Để đưa nhiều đồ thị lên cùng một cửa sổ ta dùng lệnh subpot, chính xác hơn,
subplot(m,n,i) tạo nan bân vẽ, được sắp xếp ừong mảng với m dòng và n cộ t Lệnh này
cũng bẳt plot tiếp theo phải chuyển đến hệ tọa độ thứ i (đếm 1heo số dòng). Ví dụ:
» t=(0:J:2*pỉ);
» subpìot(2,2,l)
» plot(t,sin(t))
» subpĩot(2,2,2)
» pỉot(t,cos(t))
» stibpỉot(2,2,3)
» plot(t,exp(t))
» subpỉot(2,2,4)
» pỉot(t,l./(l+ tM ))
22
4.2 Vẽ đường
Hàm plot được sử dụng để vẽ dữ liệu trên mặt phẳng. Cho vec tơ X gồm các
hoành độ x l .,,xn và vec tơ y gồm các tung độ y l,...yn, lệnhplot(x,y) sẽ vê các điểm
từ (xl,yl) đến (xn,yn). Theo mặc định, các điểm này sẽ được nối theo thứ tự bởi các
đoạn thẳng.
Câu lệnh tổng quát:
Mâu của dường vẽ Ký tụ hiên thị Nét vẽ

y yellow . dâu châm(point) solid
m magenta 0 vòng tròn : dotted
c cyan X dấu X dashdot
r red + dấu cộng dashed
g green * dấu sao
b blue s hĩnh vuông
w white d kim cương
k black A tam giác chi lên
V tam giác chỉ xuống
< tam giác chi sang
trái
> tam giác chi sang
phải
p ông sao 5 cánh
h ông sao 6 cánh
Ví dụ:
Lệnh plot(x,y, ’ro: ’) vê đường nối bởi các chấm màu đỏ với điểm dữ liệu là vòng
tròn nhỏ.
Lệnh plot(x,y,’y-’,x,y,’go’) vẽ dữ liệu 2 lần với đường nối liền nét màu vàng và
các điểm dữ liệu là các vòng tròn xanh lá cây.
Để vê tứ giác với các đỉnh (0,0), (1,1), (4,2), và (5,-1) ta có thể dùng lệnh:
» x=[0 1 4 5 0];
» y = f 0 1 2 - 1 0];
» plot(x,y)
Ta cũng có thể gần đúng đồ thị hàm số bằng cách lấy nhiều điểm. Ví dụ vẽ đồ thị
hàm số y = x 3/ [-2;2]. Trước hết, ta lấy vec tơ dòng với các thành phần biến thiên từ -2
đến 2 với độ dài bước chia là 0.05
» x=-2:.05:2;
Tiếp đó là xảc định các giả trị của y bằng cách lấy lũy thừa tùng phần tử của X.
» y=x.A3;
Cuối cùng:
» Plot(x,y)
Ta đặt tên cho hình vẽ nhờ đánh lệnh
» titỉe(‘D o thi cua ham f(x)=xA3 ’)
Do thi cua ham so ((?:)=X3
Hình J_i ĐỒ thị hàm số f i x ) = X3

Vẽ cảc đường cong tham số cũng tương tự:
r(t) = (2í cosí/(í + l),2í sừư / ự +1)) với t e [0,4#],
» t=0:J:4*pi;
» x=2*t. *cos(t)./(t+ỉ);
» y =2 *t. *sin(t)./(t+l);
» pỉot(x,y);
» titỉe( '(2t COS i/ (t+l)),2tsin t/(t+l)f)
Chủ v: các phép toán nhân chia phải được ihục hiện trên từng phần tử(nghĩa là
./ hoặc .A).
Để có tỷ lệ chính xác ta dùng:
» axis equal
24
05 Iỉ
Hình 2 : Đồ thị hàm sổ cho bởi phương trình tham số
Để vẽ nhiều đường trên cùng 1 hình ta dùng lệnh “hold ori\
Ví dụ: vẽ 2 đường tròn X1 + y 1 = 4 và (x -1 )2 + ( y - l ) 1 = 1. Chuyển sang phương
trình tham sế:
rxự) = (2cos t,2 sin í) và r2(í) = (1+cos t,1+ sin í) với t € [0,2#].
» t=ồ:pi/20:2 *pi;
» pỉot(2*cos(t),2*sin(t))
» hoỉd on
» pỉot(ỉ +cos(t),l +sin(t))
» axỉs equaỉ
» titỉeCCac dtiong tron xA2+yA2=4 va (x-ỉ)A2+(y-ỉ)/x2 = ĩ')
26 -2 15 -1 -05 ữ 0,5 1 5í ỉ 25
Hình 3: vẽ đồ thị nhiều hàm trên cùng một hình

Trong không gian 3 chiều thì thay lệnh plotbằng plot3. Ví dụ:
Để vẽ đường tham số r(í) = (cos(0»sin(í),í) với t € [0,8^] ta có thể làm như sau:
t=0:.ỉ:8*pi;
» pỉot3(cos(t),sin(t),t)
25
» titỉe(Ỵcỡs t,sin ự ) ’)
(coĩt,sin1,1}
Hình 4a: vẽ đồ thị đường cho bởi phương trình tham số trong không gian 3 chiều
Đê vẽ đường trong tọa độ cực p = 0 2, vởi 0 < 0 < n :
» iheta = 0:.2:5*pi;
» rho=theta.A2;
» poỉar(theta,rho,fhb')
4.3 Vẽ mặt
Để đồ thị hàm số / (x,y) trên miền hình chữ nhật
R =[a,5]x[c,í/] = {(xSiy)|a< x<b& c < y < d },
Trước hết tã tạo lưới điểm trong miền khảo sát nhở hàm meshgrỉd
Ví dụ: vẽ đồ thị hàm số f(x, y) = x2y - 2y ừên miềa chữ nhật [-2,2] X [-1,1]. Ta sử
dụng lưới vuông vói độ dài bước chia 0.1
>> [x,y]=meshgrid(-21: 2 , - ỉ ỉ :1);
Sử dụng phép toán vec tơ hóa để xác định z.
26
> > z=(x. A2). *y-2 *y;

Cuối cừng ta vẽ
» surf(x,y,z)
» ũtỉeựĐo thi cua ham so f(x,y) =xA2y-2y t)
Do thi cua ham »0 f(*,vì=»"y-2y
Hình 5: Vẽ mặt f ( x ty) = x 2y - 2y ữong 3 chiều

Một trong những khó khăn khi vẽ mặt là phải đối mặt với phép chia cho 0. Chẳng
hạn, ta muốn vẽ đồ thị của hàm số
/ (x, y) = . ^ trên lưới vuông [-1,1] X[-1,1] • Thực hiện lệnh

V*2+y2
>> ỉ'x,y]=meshgrid(-l:J:l);
» z=x. *ỵ./sqrt(x.A2+y.A2);
Nếu ta cứ vẽ mặt như vậy điểm lưới tại (0,0) sẽ khuyết
Hình 6a: Vẽ măt fipc, y) = . xy
27
Đê tránh khỏi tình huống này, một cách đơn giản là hãy xác định lưới sao cho
điểm (0,0) không thuộc điểm chia. Chẳng hạn:
» surf(x,y,z)
» [x,y]=meshgrid(-.99:.l:l);
» z=x. *y./sqrt(x. A2 +y. A2);
» surf(x,y,z)
» titỉeựDo thi ham so f(x,y)=xy/\surd(xA2Jry^2ỵ)
» axis equaỉ
DoIhi hamso
Hình 6b: Vẽ mặt ừong 3 chiều tránh khuyết

Ví dụ. Vẽ z=sin(r)/r với r = ^ x 2 + y1 ,- 8 < X< 8, - 8 < ;; < 8.
» x=-8:.5:8;y=-8:.5:8;
» [x,y]=meshgrid(x,y);
» r=sqrt(x.A2+y.A2)+eps;
» z=sin(r)./r;
» mesh(x,y,z)
28
Điểm khác nhau giữa hai lệnh surf và mesh là: surf tô màu bề mặt,còn mesh thì
không.
29
Chương 5: Các ỉuồng điều khiển
5.1 Các toán tử logic
Lệnh Kêt quả

A=(b>c) A=1 nêu b>c, tương tự < <=
A=(b=c) A=1 nêu b “ C
A=(b~=c) A=1 nêu b khác c

<
A=1 nêu b“ 0
II
Ị
A=(b&c) A=1 nêu b =true và c=true

A-(b|c) A=1 nêu b—true hoặc c—true
Bans 1: mối quan hệ giữa các phép toán logic
Chú ý : toán tử & và toán tử I có mức ưu tiên như nhau ừong Matlab,điều đó có
nghĩa là nó sẽ thực hiện từ trái qua phải.
5.2 Lệnh find

Chúng ta có thể trích các phần tử ừong vec tơ và ma trận thỏa mãn những điều
kiện nào đó, bằng toán tử xác định vị trí. Chẳng hạn như dùng lệnh xác định vị trí find
» j c = /7 1 3 4 1 ] ;
» i = (x == 1)
i -
1 1 0 0 1
» y = x(i)
y =
1 1 1
» j =find(x = = 1) % ì là chỉ số của phần tử thỏa mãn X = = 1
j =
1 2 5
» z = X0)
z=
1 1 1
Một ví dụ khác :
> > X = -1:0.05:1;
» y = sin(x) .* sin(3*pi*x);
30
>> plot (x,y, ’) ; hold on

» k - fin d (y <= -0.1)
9 10 11 12 132930313233
» plot (x(k), y(k), ’ro ’);
» r =find (x > 0.5 & y > 0)
35 36 37 38 39 40 41
» plot (x(r), y(r), V* ’);
Dùng toán tử find trong ma trận một cách tương tự:
»A= [l 3 -3 -5; -1 2 -1 0; 3 -7 2 7];
» k =fin d (A >= 2.5)
» A(k)
arts =
Chúng ta cũng có thể xác định vị trí theo hàng cột.
» [I.J] =fin d (A>= 2.5)
» [AỢ(1), J(l)), A(I(2), J(2)), Ả(I(3), J(3))] % danh sách các giá ừị
ans =
31
Ky thuat tinh toin trong cong nghp h6a hpc
3 3 7
5.3 Cku If nh if
* )If ....E n d .
-------- ¿>ynEias --------- ----------- ¿.xcunpie -----------
i f l o g i c a l .e x p r e s s i c n i f (a > 0)
s ta te m e n t 1 b = a;
s ta te m e n t2 d is p ( ’ a i s p o s i t i v e 0 ;
.... end
end
*)If... else ...end

-------- Syntax -----
_______________ Exam ple __________
i f lo g ic a l .e x p re s s io n
block of sta te m e n ts i f (tem p era tu re > 100)
evalu ated , i f TRUE d is p ( ’Above b o i l i n g ;
e ls e to o h ig h = 1;
block of 3-t^taments e ls e
e v a lu a te d i f FALSE d is p ('T em p eratu re i s OK.’ )
end toohijgh = 0;
end
*) I f ... elseif... else ... end.
_____________ Syntax __________ ______ Exaunple ____

i f lo g lc a l_ e x p r e s s io n 1 i f ( h e ig h t > 190)
b lo ck o f s ta te m e n ts e v a lu a te d d is p (*v e ry ■ call’ ) ;
i f l o g ic a l_ e x p r e s s io n i i s TRUE e l s e i f (h e ig h t > 170)
e l s e i f lo g ic a l_ e x p r e s s io n 2 d is p ( ’t a l l ?) ;
b lo c k o f s ta te m e n ts e v a lu a te d e l s e i f (h e ig h t < 1B0)
i f lo g ic a l.e x p r e s s io m 2 i s TRUE d is p (-'sm a ll-');
e ls e e ls e
b lo c k o f s ta te m e n ts e v a lu a te d d is p ( ’ a v e ra g e 1) ;
i f no o th e r e x p re s s io n i s TRUE end
end
Chu y: bM tip dudi day c6 ihk sit dyng m-iile Script B<3 \k tap tin dupe luu b
ngoai. Chung ta c6 th£ su dyng b it klchuong trinh soan thao nao nhung phailuu vdi
duoi m o rpng la “ m” (se dupe tim hieu kl trong chuong sau).
BM tip
H Sy cho b i l t k i t qu& c u a : m , t, h , y
1. i f n > 1 a)n = 7m = ?
m = n + 2 b)n = 0 m = ?
32
else c) n = -7 m = ?
m = n-2
end
2. if s <= 1 a) s = 1 1 = ?
t = 2zb)s = 7 t = ?
elseif s < 10 c) s = 5 7 1 = ?
t = 9 - z d ) s = 300 t = ?
elseif s < 100
t = sqrt(s)
else
t =s
end
3. if t >= 24 a) t = 50 h = ?
z = 3t + l b ) t = 1 9 h = ?
elseif t< 9 c) t - -6 h = ?
z = tA2 / 3 - 2 t d ) t = 0 h = ?
else
z —-t
end
4 . i f 0 < x < 7 a) X = -1 y = ?
y = 4xb)x = 5y = ?
elseif 7 < x < 5 5 c ) x = 30 y = ?
y = -10xd)x = 56y = ?
else
y = 333
end
Một loại cấu trúc điều kiện khác là switch. Câu lệnh này cho phép thực hiện rẽ
nhánh dựa trên các giá trị biểu thức.
33
_______________ S y n t a x ____________ _______________ E xa m p le ____________

switch expression method = 2;
case choicei svitch method
block of conmandsl case -[1,2}
case -{ch.oice2a. choice2b... dispC’ Method ia linear.’ );
block of conmands2 case 3
dispC’ Method is cubic.1);
otherwise case 4
block of commands disp(’ Method is nearest.’ );
end otherwise
dispC’ Unknown m e t h o d . ’
);
end
Cấu trúc này sê ừánh sự dài dòng và khó hiểu của cấu trúc if...end khi sử dụng
quá nhiều điều kiện.
5.4 Cấu trúc lặp

*) vòng lặp for
S y n ta x Example
f o r in d e x “ f i r s t : s t e p : l a s t s u m x = 0;
b lo ck of B ta te n e n ts for i ” l :le n g th (x )
end Bumx = sum x + x ( i )
end
Chỉ Số của phép lặp có thể là vec tơ:

___ C ix a m p íe J. ___ ____ c x a m y i e <5 ____ ____ C I a i r p i s đ ____ r, T.íimp I.p 4
for i-t:2:n for i-n:-l:3 for 1-0:0.6:4 for X - [26 9 81]

dísp(x"2); disp(sqrt(x));
end end end end
*) Vòng lặp while

Các câu lệnh trong thân của vòng lặp while sẽ được lặp lại chừng nào biểu thức
expression còn là true:
--------- Dĩ/ncax ---------- ------------ £ ,xa m p i,e -------------
w h ile e x p re ssio n N = 100;
sta te m e n t1 i t e r = 1;
sta te m e u t2 msum 0;
S ta te m en ts v h i l s i t e r <= K
¿ s u m “ msuni + i r e r ;
end ite r = ite r + 1;
end:
34
Chương 6: Phương pháp sổ
6.1 Đường cong hồi quy

Vói lệnh polyíĩt bất kỳ đa thức nào có thể khớp nhau đối với dữ liệu, polyfit ( X,
y, n) để tim thấy các hệ số của một đa thức với bậc n mà phù họp với dữ liệu nhất(tìm
mối quan hệ tuyến tính giữa X và y). Chứng ta hãy bắt đầu làm việc với hồi quy tuyến
tính của dữ liệu nào đó:
>> X = 0:10;
» y = [-.10 .24 1.02 1.58 2.84 2.762.99 4.05 4.83 5.22 7.51]
» p = polyfit (x, y, 1) % tìm đa thức hồi quy bậc 1
p = 0 .6 7 7 2 - 0 .3 9 1 4
Biến ra của lệnh polyfit là 1 vec tơ hàng các hệ số của đa thức hồi quy. Như ở ví
dụ trên, hàm hồi quy tìm được là: y = 0.6772x - 0.3914.
6.2 Phép nội suy

Cách đơn giản nhất để khảo sát một phép nội suy là vẽ đồ thị hàm số bằng hàm
plot. Các số liệu lân cận được nối với nhau bởi đoạn thẳng:
>> x l = ỉinspace(0,2*pi,2);
» x2 = linspace(0,2*pi,4);
» X3 ~ linspace(0,2*pi,16);
» X4 = ỉinspace(0,2*pị,256);
» pỉot(xl,sin(xl),x2,sin(x2),x3,sin(x3),x4,sin(x4))
» legend(’2 points ‘4 points ’16points ’256points ’)
Trong Matlab cũng có những hàm nội suy như: iníerpl(hồi quy 1 chiều), interp2,
spline. Ta thực hiện lệnh sau:
» N = 50;
» X = linspace(0,5,N);
» y = sin(x) . *sin(6*x);
» subplot(2,l,l);pỉot(x,y);
» hold on
» p = randperm(N);
» pp = p ( l :round(N/2)); % chọn ngẫu nhiên N/2 số ừong [0,5]
» pp —sort(pp); % phân loại phần tử
» XX = x ( p p ); % c h ọ n đ iể m
» y y - y(pp);
» pỉot(xx,yy, ’ro-’)
35
» yn = interpl(xx,yy,x, ’nearest’) ;
» plot(x,yn, ’g ’)
» axis tight
» legendựOriginal’, ’Crude version ’Nearest neighbor interpolation ’)
» subplot(2,l,2);plot(xx,yy, ’ro-’) ;
» hold on
» yc = interpl(xx,yy,x, ’linear’) ;
» plot(x,yc, ’g ’)
» ys = spline(xx,yy,x); % hồi quy spline
» plot(x,ys, ‘k ’)
» axis tight
» legendựCrude version ’Linear interpolation 'Spline interpolation ’)
» [X,Y,Z]= peaks (10);
» [Xi.Yi] —meshgrid(-3:.25:3,-3:.25:3);
» Zi - interp2(X,Y,Z,Xi,Yi);
» mesh(Xi,Yi,Zi);
6.3 Giá trị của hàm sỗ

Các hàm khảo sát hàm số
Lệnh Kêt quả
Fplot(‘f ,[mm_x Vẽ đô thị hàm sô trên đoạn [mm_x max_x]
ma xx ])
Fmm(‘f ,min_x, Cho giá trị của X khi hàm sô đạt cực tiêu /[min x,
maxx) maxx]
Fmax(T, mm_x, Cho giá trị của X khi hàm sô cực đại /[mm_x,
maxx) ma xx ]
Ví dụ:
» x = 0:pỉ/8:2*pi;
» y = sin(8*x);
» plot(x,y, ’b ')
» hold on
» jp lo t(’sin(8*z) ’,[0 2*pi], ‘r ’)
» title(’sin(8*z)’)
36
>> hold o ff
Chúng ta cũng cỏ thể định nghĩa hàm
> > / = ’sỉn(8*z)
» jplot(f,[0 2*piJ, ’r ‘)
» title(j)
6.3.1 Hàm 'inline’

Nó có thể rất hữu ích khi ta định nghĩa 1 hàm ừong suốt thời gian chạy Matlab.
> > / = inline( ’cos(x). *sin(2*x) ’)
f =
Inline function:
f(x) = cos(x).*sin(2*x)
» g = inline(’sqrt(x.A2+y.A2) ’x \ ’y ’)
g =
Inline function:
g(x,y) = sqrt(x.A2+y.^2)
Chúng ta có thể xác định giá trị của hàm số một cách thông thường:
>>f(-2)
ans =
-0.3149
» g(3>4)
ans
5
»A= [l 2; 3 4];
» B = [2 3; 4 5];
» g(A,B) % hàm cũng cổ thể thực hiện với mảng
ans
2.2361 3.6056
5.0000 6.4031
1 1
Chăng hạn xét hàm sô f ( x) =-------- —:------------------------------------------------- + -^ , điêm cực ti
(x —0.1) +0.1 (je —1) +0.1
hàm f có thể tìm nhờ hàm ‘fmỉn’. Theo mặc định sai số tương đối là Ỉ 0 ^ , tuy nhiên, ta
có thể thực hiện với độ chính xác hơn.
>> format long % thay đổi định dạng để tăng thêm độ chỉnh xác
» / = inline ( ’l./((x-0.1).A2 + 0.1) + l./((x-l)M + 0.1)’) ;
» jp lo t( f,[ 0 2]);
» xm l - fmỉn(f,0.3,1); fin í =f(xml);
» xm2 =finỉn(f,0.3,1,[0,le-8]); fin2 =f(xm2);
» [xml, xm2J % so sánh 2 kết quả
6.4 Phép tính tích phân và tích phân

Để lấy tích phân, ta có thể dùng hàm trapz :
> > X = 0 : 0 .5 : 1 0 ; y - 0 .5 * s q r t( x ) + X . * s in (x );
» integral1 - trapz(x,y)
integral1 -
18.1655
» X = 0:0.05: ì0;y = 0.5 * sqrt(x) + X .* sin(x);
» integran = trapz(x,y)
integran =
18.3846
Có thể thu được kết quả chính xác nếu sử dụng hàm quad, hoặcquadl, vẫn lấy ví
dụ: f ( x) =-------------------------------------------------------- \ ---------- + ------------------ ----------
(jc- 0.1) +0.1 { x - \ ý +0.1
» /= '1 J((x-0.1).A2 + 0.1) + ỉ./((x-l).A2 + 0.1)’;
» integraỉl = quad(f,0,2)
» integran = quadl(f,0,2)
6.5 Tính toán sổ học và cấu trúc điều khiển

Chúng ta có thể gặp những bài toán nhu khai triển Taylor.Chúng ta có thể khai
triển hàm số — =ỉ + x + x 2 + x3 + ... + 0(jt4) , với x=0.42
l-x
s = 0; x = 1; s = 0; X = 1;
X 0 = 0 .4 2 ; xO = 0 .4 2 ;
while (x > le-6)

hoặc
while (x > (le-6)*s)
s = s + x; s = s + x;
X - X * xO; X —X * xO;
end
Sẽ cho ta kết quả gần đúng của hàm / = — .

1 -x
38
Chương 7: Viết chương trình trong Matlab
7.1 Kịch bản và hàm m-file

Kịch bản m-file
Các lệnh của Matlab được nhập sau dấu nhắc lệnh trong Matlab. Khi bài toán ừở
lên phức tạp hơn thì cách này tỏ ra không hiệu quả. Có 1 cách giải quyết, đó là sử dụng
m-file. Lưu ý là m-file phải được đặt ở một trong những thư mục mà Matlab sẽ tự
động tìm kiếm m-fĩle trong đó; danh mục các thư mục như thế có thể xem nhờ lệnh
path. Thư mục mà Matlab luôn khảo sát đó chính là ‘Current Directory’, người dùng
có thể thay đổi.
Một m-file có ửiể được biên soạn như sau:
Chẳng hạn chúng ta lưu file sinplot.m có nội dung:
X = 0 : 0 .2 : 6 ;
y - sỉn(x);
plot(x,y);
title ('Plot o fy = sin(x) ’);
và để khởi chạy nó ta dùng lệnh
>> sinpỉot
7.2 Hàm m-file

Cấu trúc
function [outputArgs] =function_name (inputArgs)
biến ra đặt trong dấu [] và biến vào đặt ừong dấu 0
Chú ý: tên của hàm và tên m-file lưu phải đồng nhất
Ví dụ: tạo hàm average có tên được lưu là average.m. Những dòng chú ý đặt sau
dòng function và trước các dòng lệnh nội dung là dòng chú thích về hàm, nó có thể
được hiển thị nhờ lệnh » help average. Có thể gọi hàm average.m bằng cách
avrl = average(l:10)
function avr = average (x)
%ÁVERAGE tỉnh toán giả trị trung bình của các giá trị trong vec tơ
% và kết quả được gán vào avr
% Chú ỷ: đây là 1 VI dụ về hàm
n = length (x);
avr = sum(x)/n;
return;
Hay 1 ví dụ khác:
39
function [avr.sd] = stat(x)

%STAT thống kê đơn giản.
% tính toán giá trị ừvng bình và độ lệch tiêu chuẩn của vec tơ X.
n = length (x);
avr = sum(x)/n;
sd = sqrt(sum((x - avr).A2)/n);
return;
7.2.1 Những biến đặc biệt của hàm số

Mỗi hàm số đều có 2 biến bên trong: số lượng đối số vào được gọi ừong hàm là
nargin và nargout là số biến ra. Việc sử dụng nargin cho phép ta đặt giá ừị ngầm
định cho các biến vào vắng mặt.
function [outl,out2J = checkarg (Ìnl,in2,in3)
%CHECKARG mô tả cách sử dụng biến nargin và nargout.
i f (nargin == 0)
ảispựno input arguments ’) ;
43
return;
elseif (nargin == 1)
s = ini;
p = ini;
disp(’l input argument’) ;
elseif (nargin — 2)
s = inl+in2;
p - ini *in2;
disp(’2 input arguments ’);
elseif (nargin - - 3)
s —inl+in2+in3;
p - inl*in2*in3;
disp(’3 input arguments ’) ;
else
error(’Too many inputs. ’) ;
end
i f (nargout == 0)
return;
eiseif (nargout == ĩ)
outl - s;
else
outl = s;
out2 - p ;
end
7.2.2 Biến địa phương và biến toàn cục

Các biến thông thường đều được mặc định là biến địa phương, nó chỉ được tham
chiếu trong hàm chứa nó. Còn các biến cục bộ được khai báo sau từ ‘global’, nó sẽ
được tham chiếu trong toàn bộ modun.
7.2.3 Tính giá trị hàm sổ cách gián tiếp

Việc tính hàm số gián tiếp giúp cho việc lập ừ ình ừở lên tổng quát hơn. Trong
Matlab ta sử dụng lệnh ‘fevaV.
[yl,..,ynj =feval (Fjclr..,xn),
Trong đó F là tên hàm số được định nghĩa , xl,...,xn là biến vào và y l,...,yn có
thể là biến ra. Xem ví dụ:
> > X = p i; y = co s(x );
» z = feval(’cos’,x);
Dòng lệnh cuối có thể được thay bởi 2 dòng lệnh
» F = ’cos’;
» z = feval(F,x)
Việc tính giá trị hàm gián tiếp sẽ là 1 công cụ tốt để xây dựng 1 chương ữình với
hàm số được coi như là tham số.
Ví dụ tạo hàm có tên là funplot
function fimplot (F, xstart, xend, col);
%FUNPLOT makes a plot o f the function F at the interval [xstart, xendj.
% The plot should be made in one o f the standard Matlab colors, so
% ’co l’ is one o f the following value: or V’.
% default values:
% [xstart,xend] = [0,10]
% co ỉ= Ẽb ‘
% Note: illustrates the use o f feval command
i f (nargm == 0)
error ( ’No function is provided. ’) ;
end
41
i f (nargin < 2)
xstart - 0;
xend= 10;
end
i f (nargin - - 2)
error ( ’Wrong number o f arguments. You should provide xstart and xend. ’) ;
end
i f (nargin < 4)
col = ’b ’;
end
i f (xstart == xend),
error ( ’The [xstart, xend] should be a non-zero range. ’);
elseif (xstart > xend),
exchange = xend;
xend —xstart;
xstart - exchange;
end
switch col
case { ’b ’, ’k ’, ’m ’, 'g ’w ’, ’y ’, V '}
; % do nothing; the right color choice
otherwise
error ( ’Wrong col value provided. ’)
end
x - linspace(xstart, xend);
y - feval(F,x);
plot (x,y,col);
description - [ ’Plot o f ’, FJ;
title (description);
return;
42
Chương 8: Văn bản
8.1 Chuỗi kí tự
Matlab hầu như làm việc vói số, nhưng khi cần thiết ta có thể xử lý văn bản. Văn
bản được lưu dưới dạng các kí tự chuỗi. Chuỗi là vec tơ mà phần tử là các giá trị
ASCn được hiển thị dưới dạng chuỗi kí tự. Do đó, các thao tác trên chuỗi giống như 1
vectơ
» t - ’This is a character string’
t=
This is a character string
» size(t)
ans
I 27
» whos
Name Size Bytes Class
11x27 54 char array
» t ( l 0 :1 9 )
ans =
47
character
» t([2,3,10,17])
ans
hi t
Để chuyển đổi dạng chuỗi về dạng số ta dùng double
» double(t(l:12))
ans -
84 104 105 115 32 105 115 32 97 32 99 104
Và hàm char được cung cấp để biến đổi ngược lại
» t([16:17])
ans -
ct
» t([16:17])+3 % thực hiện các thao tác toán học với mã ASCII
ans -
43
102 119
» t([16:17])-3 %gọi lại mã ASCII
ans
96 H ỉ
» char(t([16:17])-2) % chuyển mã ASCII ra thành kí tự
ans
ar
8.2 Vào - ra dữ liệu

Các lệnh liên quan đến nhập dữ liệu từ bàn phún: input, keyboard, menu, pause
8.2.1 Vào ra dữ liệu từ bàn phím

8.2.1.1 Lênh ‘input’
Lệnh này đưa ra thông báo nhắc nhở người sử dụng nhập dữ liệu từ bàn phún.
Lệnh có dạng: r= input(string)
Trong đó r là tên biến, string là sâu ký tự chứa thông báo nhắc người dừng biết về
dữ liệu cần nạp. Dữ liệu đó có thể là biểu thức bất kỳ. Nếu người sử dụng ấn ‘Enter’
ngay thì r sê là 1 ma ừận rỗng.
Ví dụ:
>> m=input('nhap so dong cua ma ừan: ’)
nhap so dong cua ma tran: 3
m
3
» n-ìnputựnhap so cot cua ma tran: ')
nhap so cot cua ma tran: 5
n
5
8.2.2 Đưa dữ liệu ra màn hình

8.2.2.1 lênh 'disv’
Lệnh disp để hiển thị nội dung cùa 1 mảng ra màn hình mà không hiển thị tên
mảng. Lệnh này luôn thực hiện cho dù có dấu ở cuối dòng lệnh hay không, ngoại
trừ trường hợp mảng rỗng. Lệnh có dạng:
Disp(x)
Trong đó X có thể là tên biến mảng hoặc biểu thức bất kì.
>> disp(’This is a statement. ’) % 1 chuỗi
This is a statement.
» disp(rand(3)) %1 ma trận
44
0.2221 0.0129 0.8519

0.4885 0.0538 0.5039
0.2290 0.3949 0.4239
8.2.2.2 lênh ‘Ịprintí‘
Lệnh có cấu trúc
Fprintf(dmlog_format, danh sách biến)
Trong đó dialog_format là biến(hằng) xâu ký tụ chứa thông báo và các ký tụ định
khuôn dạng dữ liệu ra.
Ví dụ:
>> a=rand(3,3)
a
0.8462 0.6721 0.6813
0.5252 0.8381 0.3795
0.2026 0.0196 0.8318
>> fprintf('Do dai cua ma tran A: %i%i\ỉength(a))
Do dai cua ma ừan A: 3
» fprintf('A bỉnhphuong: \n %i%i%i\n%i%i%i\n%i%i%i\n\aA2)
A binh phuong:
1.207124e+0009.614361e-0013.503595e-001
1.145488e+000L 062870e+0001.690032e-001
1.398258e+0009.914752e-0018.373965e-001
Hay một ví dụ khác:
>> sfr = ’life is beautiýul’;
» Jprintf(’My sentence is: %s\n ’,str); % chủ ỷ định dạng \n
My sentence is: ỉife is beautiýuỉ
» fprintf(‘My sentence is: %30s\n ’,sư);
My sentence is: ỉife is beautiýuỉ
» fprintf(’My sentence is: %30.10s\n\str);
My sentence is: ỉife is be
» fprintf(’My sentence is: %-20.10s\n \str);
My sentence is: ỉife is be
»
» name - ’John
» age = 30;
45
50
» salary = 6130.50;
» Jprintf(’My name is %4s. I am %2d. My salary is f %7.2f.\n’,name, age,
salary);
» My name is John. I am 30. My salary is f 6130.50.
»
> > x = [0,0.5, 1];
» y - [x; exp(x)];
» jprintf(’%6.2f%12.8j\n \y);
0.00 1.00000000
0.50 1.64872127
1.00 2.71828183
»
» jprintf(’%6.1e %12.4e\n’,y);
0.0e+00 1.0000e+00
5.0e-01 1.6487e+00
L0e+00 2.7183e+00
»
» x = 1:3:7;
» y = [x; sin(x)J;
» jprintf(’%2d %10.4g\n \y);
1 0.8415
4 -0.7568
7 0.657
8.2.3 Xử lífile văn bản

8.2.3.1 Lênh ‘fopen’
fllelD =fopen(fĩlename, permission)
Lệnh được sử dụng để mở một file văn bản có tên được cho bởi filename ở chế độ
được chỉ ra bởi permission. Permission có thể là:
‘rt’ đọc
‘w t’ ghi (sẽ ghi 1 file mới nếu không tìm thấy filevới tên chỉ ra bởi filename)
‘at’ nối đuôi(sẽ ghi 1 file mới nếu khôngtim thấy file với tên chỉ ra bởi
filename)
‘rt+’ đọc và ghi (không tạo file)
Biến nhận dạng fîlelD sẽ nhận các giá trị nguyên và ta sẽ thâm nhập vào file sử
dụng biến nhận dạng này.
Khi kiểm tra lỗi mở file có thể sử dụng
[fid, message] =/openựUename, permission)
8.2.3.2 Lênh đóng file ‘fclose’
S t =fclose(fîd)
Lệnh này sẽ đóng file ứng với biến nhận dạng fid. Fclose sẽ ừả biến st giá trị 0
nếu nó là thành việc đóng file và trả lại giá trị -1 nếu gặp lỗi.
Lệnh st = fclose(‘aU’) đóng tất cả các file đang mở ngoại trừ 0, 1 và 2.
8.2.3.3 Lênh ‘frewind’
Frewind(fid)
Đặt con ừỏ file có tên fid vào đầu file.
8.2.3.4 Lênh ‘fscanf’
[A, count] =fscan/(fĩd, format, size)
Lệnh này đọc dữ liệu từ file tương ứng với biến nhận dạng fid, chuyển đổi dữ liệu
về dạng được xác định bởi xâu biến format, và gán vào ma trận A. Biến count là biến
ra tùy chọn dùng để chứa số lượng phần tử đọc được.
Fid là biến nhận dạng tương ứng với file nhận dạng được từ lệnh fopen.
Size là tùy chọn; giới hạn số đối với số phần tà được đọc tù file, nếu không có
biến này thì toàn bộ file sẽ được xét Các giá trị của biến có thể là:
• N: đọc không quá N phần tử từ file vào vec tơ cột
• Inf: đọc không quá kết thúc file
• [M,N] : đọc không quá M*N phần tử và đưa vào ma trận kích thước
không quá M*N theo từng cột. N có thể là inf còn M thì phải là hữu hạn
Nếu ma trận A là kết quả của việc chuyển định dạng kí tự và biến size không có
dạng [M,N] thì vec to dòng sẽ được trả lại.
Format là biến xâu chứa các kí tự chuyển đổi định dạng của ngôn ngữ c . Các kí
tự định khuôn dạng bao gồm : %, các kí tự thay thế, độ dài trường, và các kí tự chuyển
đổi dạng: d, i, o, u, X, e, Ẹ, g, s, c and [... ](liệt kê tập hợp).
Nếu %s được sử dụng thì khi đọc 1 phần tử có thể dẫn tới 1 loạt các thành phần
của ma trận được sử dụng, mỗi thành phần giữ 1 kí tự.
Sử dụng %c để đọc kí tự trắng; định dạng %s bỏ qua các kí tự trắng.
Neu chỉ thị định dạng gồm lẫn lộn cả số lẫn kí tự thì ma trận kết quả sẽ là ma trận
số và mỗi kí tụ sê chuyển thành 1 số chính bằng giá trị mã ASCII của nó.
Fscanf khác với lệnh này trong ngôn ngữ c ở chỗ nó là lệnh được vec tơ hóa để
ừả lại đối số là ma trận. Biến xâu định dạng sẽ được lặp lại cho đến khi gặp kết thúc
file hoặc đọc đủ số lượng phần tù chỉ ra bởi size.
Ví dụ:
47
Lệnh
s =fscanf(fỉd,’%s’)
Đọc(và ừả lại) một xâu
Lệnh
Á =fscanf(fỉd,’%5d’)
Đọc các số có 5 chữ số thập.
Ví dụ: giả sử có file văn bản với tên ‘kq.txt’ chứa xâu ‘Day la ket qua dua ra’. Khi
đó ta có thể đọc dữ liệu vào như sau:
>> fid - f o p e n ( ' k q ’rt’)
Fid =
3
» x= fscanf(fìd, ’%s ’)
x=
Daylaketquaduara
» frewind(fid)
»x= fscanf(fid, ’%s%c ’)
x=
Day ỉa ket qua dua ra
Yí dụ: giả sử fĩle văn bản Matrix.txt chứa 2 dòng
1 2 3 4 5
6 7 8 9 10
Hãy theo dõi sự tác động của lệnh fscanf
>>fopen(‘m a ừ ix .tx t'r t’)
Ans =
4
» A=fscanf(4, ’% i’,[2,5])
A =
12 3 4 5
67 8 9 10
>>frewind(4);
» B =fsccmf(4, ’Voi’,[5,2])
B=
1 6
2 7
48
3 8
4 9
5 10
» Frewind(4');
» C=fscanf(4, ’% i’,6)
c=
1
2
3
4
5
6
8.2.3.5 Lênh FPRINTF

Count = fprintf(fidformat,A,...)
Lệnh này sê định dạng các thành phần của ma trận A( và các biến tiếp theo trong
danh sách) theo định dạng được xác định bởi xâu format, và ghi ra file tương ứng với
các biến nhận dạng fid. Biến count sẽ đếm số byte của dữ liệu được ghi ra. Fid là số
nguyên nhận dạng tên file thu được từ lệnh fopen. Có thể dùng số 1 nếu sử dụng thiết
bị ra chuẩn(màn hình). Nếu trong lệnh này không chỉ ra biến fid thì kết quả sẽ đưa ra
màn hình.
Format là biến xâu chứa các kí tự mô tả đinh dạng dữ liệu giống như trong ngôn
ngữ c. Các kí tự \n, \r, \t, \b, \f có thể dùng để tạo ra linefeed, carriage return, tab,
backspace, và formfeed characters tương ứng. Sử dụng \\ để tạo dấu \ và %% để tạo ra
kí tự %.
Ví dụ:
>> X=0:.1:1; y=[x,exp(x)];
» ýìd~fopen(’exp.txt', 'w');
» fprintf(fid, '%6.2f %12.8J\r’,y);
» fclose(fid);
Đe tạo ra file văn bản có tên là ‘exp.txt’ chứa bảng giá trị hàm mũ:
0.00 0.10000000
0.20 0.30000000
0.40 0.50000000
0.60 0.70000000
0.80 0.90000000
49
1.00 1.00000000
1.11 1.22140276
1.35 1.49182470
1.65 1.82211880
2.01 2.22554093
2.46 2.71828183
In nội dung fĩle đó ra màn hình:

> > jprintf('**** Day la bang gia trì ham exp ****(« X exp(x)\n’) ;
**** Day la bang gia tri ham exp ****
X exp(x)
» jprintf('%6.2f %12.8f\n',y);
0.00 0.10000000
0.20 0.30000000
0.40 0.50000000
0.60 0.70000000
0.80 0.90000000
1.00 1.00000000
1.11 1.22140276
1.35 1.49182470
1.65 1.82211880
2.01 2.22554093
2.46 2.71828183
50
Phần II: ứng dụng của Matlab

trong công nghệ hóa học
Trong thực tế, các bài toán công nghệ hóa học rất phức tạp, thiết kế và vận hành
một quá trình thì luôn luôn đi kèm với tính kinh tế của nó. Do đó, tối ưu hóa quá trình
thiết kế và vận hành là vô cùng cần thiết, việc mô phỏng và tính toán nhanh và chính
xác luôn được đánh giá rất cao. Hiện nay, chiếc máy tính cùng với các phần mền tính
toán và mô phỏng giúp chúng ta giải quyết những bài toán đó đơn giản hơn, thời gian
được tiết kiệm triệt để.
Trong học tập, các bài toán công nghệ hóa cũng cần một khối lượng tính toán khá
lớn, có rất nhiều hằng số phải tra cứu(khối lượng riêng, nhiệt hóa hơi, hằng số ừong
phương trình Antoine...) và cũng không ít những phương ừình thực nghiệm phức tạp
ừong các trường hợp khác nhau( tính chuẩn Nuselt...), việc tính toán bằng tay ừở nên
mất thời gian và có khả năng sai rất lớn. Matlab bằng công cụ tính toán, và mô phỏng
phong phú và rất mạnh có thể khắc phục điều đó.
51
Chương 1: Dẩn nhiệt và đếỉ lưu
1.1. Dần nhiệt
1.1.1.Tóm tắt li thuyết

dẫn nhiệt qua vách phẳng
ỲỂl w /m 2
1Ậ
q mật độ dòng nhiệt q = const qua các lớp 5; chiều dày của các lớp thứ i (m);A,i là
hệ số dẫn nhiệt của lớp thứ i (W/mK); tw! nhiệt độ mặt ừ o n g ; t^n+i) nhiệt độ mặt ngoài
của lóp thứ n
dẫn nhiệt qua vách trụ
tK x ^w(n+l)
«/= —
ln ^
! lĩú ị dt
trong đó :
+ qi mật độ dòng nhiệt trên một mét chiều dài (W/m)
+ di đường kính của lớp thứ i (m)
+ hệ số dẫn nhiệt của lớp thứ i (W/mK)
dẫn nhiệt không ổn định
với tâm phang dày 20
nhiệt độ tại tâm tấm
0 « =f(Bi/Fo)
Nhiệt độ trên bề mặt tấm
ex_, = f(Bi/Fo)
Trong đó :
jy. aö
Ol = — tiêu chuân Biot
ar , _
g ĩ tiêu chuẩn fourier ( a ~ £ hệ số nhiệt độ )
v _ X
$'• kích thước không thứ nguyên
52
1.1.2. Vi dụ
1. Ví dụ vách phẳng hai lớp có bề dày và hệ số dẫn nhiệt tương ứng là ỗi = 10 cm
Ằ]= 2.5W/m độ ô2=0.3m Ầ2=1.5W/m độ. nhiệt độ mặt phải là 1^2=25° c khi có dòng
nhiệt q = 500W/m2 dẫn qua vách xác định
a)nhiệt độ mặt ừái tjni và nhiệt độ mặt tiếp xúc
b) gradiel tại mỗi lóp
c)nếu giữ nguyên lớp có gradien nhỏ và duy ừì gradt như cũ thì lớp còn lại phải
thay đổi độ dãy V và 5’bằng bao nhiêu để gradt nhu nhau trên cả hai vách khi nhiệt độ
các mặt và dòng điện không đổi
Lời siải
a) xác định nhiệt trở của mặt trái t2
nhiệt trở của các lóp :
lẠpl: J = | ỉ = 0-04" 2đb/ ^ 1 ớ P 2 : ^ = y | = 0-2«12^ / w
nhiệt trở tổng R = 0.04+0.2=0.24m2độ/w

q 500 n
đô chênh lêch nhiêt đô hai măt là At = ti - Í2 = D = 7 = 2 0 8 3 . 3 c
K U.24
nhiệt độ mặt trái là: tj = 2083.3+25=2108.3

nhiệt độ chỗ tiếp xúc : tfc = ti - qRi= 2108.3 - 500.0.04 = 2088.3°c
b) tính gradiel tại các lớp
l ớ p h g r a d t ^ J = f f = 200aCI M
lóp 2 : gradt,= | - = f y = 333.33°C/M
Giải bằns Mathlab

Nhập những dòng lệnh sau, sau đó cho chạy trên mathlab ta sẽ có kết quả
function btl
delta=mputCnhap duoi dang ma ữan gia tri cua do day (don vi m ):');
lamda = inputCnhap duoi dang ma tran gia tri cua he so dan dien (don viw /m do):');
tl=input('nhap gia tri cua nhiet do mat trong cung:');% nếu không cho thì nhập là ‘no’
t2= input('nhap gia tri cua lop ngoai cung:'); % nếu không cho thì nhập là ‘no’
q=input('nhap nhiet luong cua dong: '); % nếu không cho thì nhập là ‘no’
barrier- input('loai vach ma ba can tinh:');
R=nhiettro (barrier, lamda, delta);
if strcmp(tl,'no')
tl=t2+q/sum(R)
elseif strcmp(t2,'no')
t2=tl-q/sum(R)
elseif strcmp(q,'no')
disp(tinh q nao!')
q =density(tl,t2,R)
end
digpC......................................')
dispCgia tri cua gradiel:')
gradt=gradiel(q,lamda)
dispC...............................');
dispCnhiet do cua lop tiep xuc la :');
ttx=nhietdotiepxuc(q,tl ,R);
kết quả thu được là:

»BT1
nhap duoi dang ma tran gia trì cua do day (don vi m): [0.1 0.3]
nhap duoi đang ma tran gia trì cua he so dem đien (don viW/mdo): [2.5 1.5]
nhap gia tri cua nhiet do mat ừong cung: 'no'
nhap gia trì cua lop ngoai cung: 25
nhap nhiet ỉuong cua dong: 500
loai vach ma ba can tinh: 'ýlat'
R =
0.0400 0.2000
tl =
2.1083e+003
gia ừỉ cua gradìel:

gradt -
200.0000 333.3333
nhỉet do cua lop tiep xuc la:

ttx =
2.0883e+003
54
2. Vách trụ hai lóp đường kính tron cùng di = 20cm, bề dày và hệ số dẫn nhiệt
hai lóp tương ứng là
5i = 2cm Ằ] = 1.2 w /m độ ,52 = 3cm k2 = 0.8 w /m độ. Nhiệt độ mặt trong cùng và
ngoài cùng là tmi = 80°c, tm2 = 20°c. xác định
a)dòng điện dài qL qua vách nhiệt độ tại chỗ tiếp xúc
b) mật độ dòng tại chỗ tiếp xúc
c)Gradt tại mặt trong cùng
Lời siải
a)dòng điện dài qL qua vách nhiệt độ tại lóp tiếp xúc ttx
đường kính các lớp
d2 = di + 2Ỗ! = 0.2 + 2*0.02 = 0.24 m
d3= d2+ 252= 0.24 +2*0.03 = 0.3 m
nhiệt trở dẫn nhiệt tại mỗi lóp của vách trụ
1ỚP2 ¿ ^ ¿ ^ ( ẩ r 0-0444 ™ ^
nhiệt trở dẫn nhiệt tổng là R = Rt! + Rt2 = 0.0685 m độ/ w
At 80-20
mât đô dòng điên dài 9l =~z = n _ = 875.91 w /m
K U.Uoo5
nhiệt độ tại lớp tiếp xúc : t*i = tn,! - qLR] = 80 - 875.91*0.024 = 58.97°c
b) mật độ dòng tại lớp tiếp xúc chỗ có đường kính d2
qL 875.91
q = - ^ ~ = — — = 1 1 6 1 v ự /™ 2
*d2 3.14*0.24 w/m
c) gradt tại mặt trong cùng có đường kính di :
q _ q L 1 875.91 _ 1 1 £ 1 0
Ig r ad ti l- ^ ^ ^ 1.2^0.02
Giải bằne Mathlab

Với bài tập này ta cũng cho chạy chương ứình như trên và thu được kết quả như
sau chú ý thay đổi khi nó hỏi vách loại gì lúc này bạn phải nhập là ‘pier’( trụ)
kết quả thu được là :
»BT1
nhap duoi dang ma tran gia trì cua do day (don vi m): [0.02 0.03]
nhap duoi dang ma tran gia tri cua he so dan diert (don viW/mdo): [1.2 0.8]
nhap gia trì cua nhiet do mat ừong cung: 80
nhap gia trì cua lop ngoai cung: 20
55
nhap nhiet ỉuong cua dong: 'no'

loai vach ma ba can tinh: 'pier'
duong kinh cua ỉop trong cung: 0.2
d =
0.2000 0.2400 0.3000
R =
0.0242 0.0444
tinh q nao!
Q=
874.9653
gia tri cua gradiel:

gradt =
ỉ.0e+003 *
0.7291 1.0937
nhỉet do cua ỉop tiep xuc la:

tíx
58.8423
1.2. Đối lưu
1.2.1. Tóm tắt lí thuyết

khi tính toán lượng nhiệt ừao đổi bằng đối lưu ta dùng công thức Newton :
Q = aFAti
Trong đó :
Ati = tw- ti
Q lượng nhiệt trao đổi ừong một đơn vị thời gian là một giây
F là diện tích bề mặt ừao đổi nhiệt m2
t wlà nhiệt độ trung bình của bề mặt vật rắn
ti là nhiệt độ của môi trường
a là hệ số tỏa nhiệt W/m2K
Nu = f(Re,Pr,Gr....)
Trong đó:
a.l
Nu = ^ t i ê u chuân Nussel
a>.l ,
Re = tiêu chuẩn Reynold
Pr = — tiêu chuẩn Pradtl

a
_ g-0d3.At
Gr =---- —---- tiêu chuẩn Grashoft
V
Với:
a hệ số tỏa nhiệt w/m2 K
Ằ,hệ số dẫn nhiệt W/mK
0) tốc đọ chuyển động m2/s
a hệ số dẫn nhiệt độ m2/s
g gia tốc trọng trường m/s2
Àt = tw - ti
p là hệ số giãn nở thể tích “K"1
với chất lỏng ừa bảng
với chất khí p = 1/T
1: kích thước xác định m
dối lưu tự nhiên
tỏa nhiệt đối lưu tự nhiên trong không gian vô hạn (Re<2300)
Đối với ống hoặc tấm đặt đứng khi (Gr.Pr)>l o9
Grt = 0.15(Gr;.Pr,)° ĨĨL

Pr
Đối với ống hoặc tấm đặt nằm ngang khi 103<(Gr.Pr)<l o9
PỊỊ_
Nu, =0.5(ơrí.Pr;)°
Gr,
Nhiệt độ xác định là nhiệt độ chất lỏng hay khí ti;kích thước xác định với ống
hoặc tấm đặt đứng là chiều cao h với ống nằm ngang là đường kmh với tấm nằm
ngang là chiều rộng
đối lưu cưỡng bức khi có chất lỗng chuyển động trong ống
chế độ chảy tầng
Nu, = O.lSRe,033 Prf43 Gr,01 ĨĨ L

\ GrU
Đối với không khí
Nui =
= 0.13.Rei°'33.
0.13.Rei . (Gri°'43 '£i.eR
tỏa nhiệt khi chất lỏng chảy rối

57
/ Pr,N'025
Nu, =0.021 Re“8Pr;043
Pr
Trong đó £1 là hệ số ảnh hưởng của chiều dài ống
£r hệ số ảnh hưởng của độ cong
chất lỗng chuyển động ngang qua chùm ống
n - tổng số hàng ống do hệ số tỏa nhiệt hàng ống thứ i
i tăng di tăng nhưng i > 3 (Xi = const (a3 = a4 = a 5 = ....)
a ỉ +a1+{n-2)a3
hệ SÔtỏa nhiệt của chùm ông a = ----------- ----------
° n
để xác định a3 có
dối với chùm ống song song
( Pr, ì
Nu, = 0 .2 6 R e °65 Pr,033 €m€ c
Pr
Trong đó £9 hộ số xét ảnh hưởng góc va đập 9 , nếu 9 = 90°, e* = 1

es hệ số ảnh hưởng của bước ống theo chiều sâu
, s N0..5
đối với chùm ống so le

r ^ N0-25
Pr,
Nu, = 0.41Re"6Pr°'33
vPrw
Trong đó Eẹ hệ số xét ảnh hưởng góc va đập V , *p = 90 thì eẹ=1
^ I/6
.../Xị = 0.6a3
— > 2..............£s = 1.12...................a2 = Q.laì

s2
Si là bước ngang s2là bước dọc
Trong các công thức trên thì Re = 103 - 105. kích thước xác định là đường Tcính
ngoài. Nhiêt độ xác định là nhiệt độ trang bỉnh của chất lỏng ti
tỏa nhiệt khi có biến đối pha
tỏa nhiệt khi sôi
khi sôi bọt ở áp suất p = 0.2 - 0.8 bar
a = 46-At2V w/m K
Àt tw- ts
58
tw là nhiệt độ vách đốt nóng

Ịj nhiệt độ bảo hòa ứng với áp suất sôi
p là áp suất sôi bar
tỏa nhiệt khi ngưng màng
ngưng màng trên bề mặt vách đứng hoặc ống đứng
a t =0.943
Trong đó :
g là gia tốc ừọng trường g = 9.81m/s2
À.hệ số dẫn nhiệt của chất lỏng ngưng W/mK
r nhiệt hóa hơi J/kg
p khối lượng riêng của chất lỏng ngung kg/m3
V độ nhớt động học m2/s
h chiều cao của vách hoặc ống đứng m
d là đường kính ngoài của ống m
twnhiệt độ bề mặt vách °c
ts nhiệt độ bảo hòa ứng với áp suất ngưng tụ °c
trong các công thức trên nhiệt độ xác định là tn, = 0.5(tw- ts)
1.2.2. Vỉ dụ
1. bao hơi của lò hơi đặt nằm ngang có đường kính d = 600mm. nhiệt độ mặt
ngoài lóp bảo ôn
tw= 60°c, nhiệt độ xung quanh ti = 40°c. xác định nhiệt tỏa ra từ lm 2 bề mặt của
bao hơi tới không khí xung quanh
Lời siải
từ nhiệt độ không khí ti = 40°c tra bảng ta có X = 0.0276W/mK V = 16.69.10"6
m2/s Pr != 0.696
từ nhiệt độ t„ = 60°c ta có Prw= 0.696. ta nhận thấy Pĩ! và Prwtương đương bằng
nhau nên
(Pr,/Prw)°-25 = 1
Theo tiêu chuẩn Gr: Grt = ^^7 — At

V
ở đây thì g = 9.81m/s2, p = J, - 40+273 - ° '0032’ Át = tw -ti = 20°c
_ 9.81.0,0032.0,63.20
= 4.87.10'
r ~ (16.69.10~6)2
59
Gr !.Pr! = 4.87.10 .0,699 = 3.4.10

Do 103< PrGr < 109nên
ên áp dụng công thức
Pr,
Nu, = 0.5(Gr;.Pr,)025 = 0.5*(3.4*108)025 =68
Gr,
Mà
a.l Nu.Ẵ 68*0.0276 ^
Nu = —-=>a = —— =-----—---- = 3.13,.w/m K
Ầ l 0.6
Nhiêt lượng tỏa ra trên lm 2 bề mặt bay hơi là
q = aÀt = 3.13*20 = 62.6 w/m2
sỉải bằns Mathlab
Tạo m-file với tên là “bail.m” và có nội dung như sau:
function bail
d=input('duong kinh cua ong nam ngang: ');
twînput('nhiet do ben ngoai lop bao o n :');
tl=input('nhiet do khong k h i:');
% bang cach ừa bang ta tinh duoc nhung gia tri cua lamda ,v,Pr
[p, Cp,lamda 1,a,u,v 1,Prw] “ trab angkhikho (tw);
[p,Cp,lamda2,a,u,v2,Prl]=trabangkhikho(tl);
beta_hesogiannothetich(tl);
Grl=tieuchuanGrashoft(9.81,beta,d,v2,tw,tl);
Nul=tieuchuanNusselt( 'doiluutunhienVO', Grl,Prl,Prw,0,0,0,0,0);
disp('.............................................')
disp('he so toa nhiet la :')
alpha=Nul*lamda2/d
disp('.........................................')
disp('nhiet luong toa ra tren be mat chat bay hoi la :')
q=alpha*(tw-tl)
để giải bài tập 1 thì trên màn hình “CommandWindow” gọi » bail sẽ thu được
» bail
ảuong kình cua ong nam ngang: 0.6
nhiet do ben ngoai ỉop bao on: 60
nhiet do khong khi: 40
60
he so toa nhỉet ỉa:

alpha =
3.1260
nhiet ỉuong toa ra tren be mat chat bay hoi ỉa:

q=
62.5199
2.tính hệ số tỏa nhiệt trung bình của dầu máy biến áp chảy trong ống có đường
kính d=8mm, dài lm, nhiệt độ trung bình của vách ống tw=20°c. tốc độ chảy trong ống
là 0) = 0.6m/s. biết t[ = 80°c
Lời siải:
tra bảng các thông số phụ thuộc vào nhiệt độ của dầu máy biến áp ta có tại t =
80°c
X = 0.1056W/mKv=3.66*10"6 m2/s p = 7.2*10^0K"1, Pli = 59.3, Prw = 298
R e , = ^ = Q-6*8*10;= 1 3 1 0
' Vị 3.66*10
Rei <2300 dầu chảy tầng do vậy :
r Pr, V'25
Nu, = 0.15 Re,033 Pr,043 Gr,01
Pr
v11»;
_ g M & _ 9-81 * 7.2 *10^ * (8* 1O'3f (80 - 20) _
v2 (3.66*10^)2
1. Nui = 16.3
Nu,.A, 16.3*0.1056 ,
2. g .io - ~ 215« / A
eiài bâns mathlab
function bai2
d=input('duong kinh cua ong la d = ');
tw=input('nhiet do trung binh cua vach ong tw = ');
w=input('toc do dong chay la w = ' );
tl=input('nhiet do cua chat long tl= ');
[p l ,Cpl ,lamdal,ul ,vl ,al ,betal ,Prw]=trabangdaumaybienap(tw);
[p2,Cp2,lamda2,u2,v2,a2,beta2,Prl]=trabangdaumaybienap(tl);
Rel = tieuchuanReynold(w,d,v2);
Grl=tieuchuanGrashoft(9.81 ,beta2 ,d,v2,tw,tl) ;
Nul=tieuchuanNusselt( 'doiluucuongbuc','0', Grl,Prl,Prw,Rel,0,0,0,0);
61
disp('.................................................')
dìspChe so toa nhiet trong ong la :')
alpha=Nul*lamda2/d
>> baữ
duong kinh cua ong la d = 0.008
nhiet do trung binh cua vach ong tw = 20
toc do dong chay ỉa w = 0.6
nhiet do cua chat long tỉ= 80
he so toa nhiet trong ong la:

aỉpha
215.5734
3. một chùm ống so le gồm 10 dãy ống. đường kính ngoài của ống là d = 38mm.
dòng không khí chuyển động ngang qua chùm ống có nhiệt độ trung bình tl=500° c.
tốc độ của dòng không khí là 12m/s xác định hệ số tỏa nhiệt trung bình của chùm ống
Lời ziải i
tra các thông số ở bảng 7 ta có ứng với 500° c X - 5.74*10"2 W/mK; V =
79.38*10 m2/s; Pĩ! = 0.687
m,.ã 12.*38*10 3 , nic:
tính R-S; —— ———-------- 1T~ —5745 —>
' V, 79.38 *10“
s
Pr;
Nu, = 0.41 Re"6Pr,° eẹes = 0.41Re“6Pr;053 do bỏ qua ảnh hưỏrng của Es
PrAw y
VA vPrw/
= 1
Nu, = 0.41 *574506*0.687°33 = 65.2 => 03= ^ ^ = ^ ^ ^ ^ - = 98.5 W/m2K
a, + a 2 + (n- 2)a, 0.6a, +0.7«, + (10-2)a, ,

a = —---- — ---- — = ----- 1------ -------------— = 91.6w/m K
n 10
eiải bans Mathlab

nhập câu lệnh sau và cho chạy chương ừình trong mathlab
function bai3
n = input('so ong n = ');
d = input ('duong cua moi ong d = ');đom vị là met
tl = input('nhiet do cua dong khong khi chuyen dong ngang qua chum ong tl = ');
w = input('toc do cua ong khong khi omega = ');
62
[p,Cp,lamda,a,u,v,Prl]=trabangkhỉkho(tl);% ừa bảng khí khô để tìm ra những

thông số lamda, Pri V
Rel = tieuchuanReynold(w,d,v)
Nul=tieuchuanNusselt( 'doiluucuongbucVsole', 0,Prl,Prl,Rel, 1,1,0,0);
disp('------------------------------ ');
dispChe so toa nhiet cua ong la :')
alpha = hesotoanhiet('sole',NuỤamda,d,n)
» baiĩ
so ongn = 10
duong kinh cua moi ong d - 0.038
nhiet do cua dong khong khi chuyen dong ngang qua chum ong tỉ - 500
toc do cua ong khong khi omega = 12
he so anh huong goc va dap: 1
ephi -
1
he so anh huong buoc ong: 1
es =
1
he so toa nhiet cua ong la:

alpha =
91.6522
63
Chương 2: Bức xạ nhiệt và truyền nhiệt
2.1. Trao đỗi nhiệt bức xạ giữa các vật
2.1.1. Tóm tắt lí thuyết

hai tấm phẳng đặt song song
■^12 —Qìỉ —SqdCị A. ±2

w/m
100 100
1
Độ đen qui dẫn qd J_ + J_ _ 1
Hệ số bức xạ của vật đen tuyệt đối

Co = 5.67 w /m 2 K4
hai tẩm đật song song có màng chẳn
Khi có m màng chắn giữa với độ đen em=Ei= 62 lúc này bức xạ từ tấm phẳng 1
sang 2 sẽ giảm đi (m+1 ) lần
?12
(tm+ 1)
hai vât bọc nhau
4 ' T \4
Qn —s qdC0-^1 ■JL] 2
w
100 ì 100
1
~ £qd ~ 1 ÍT
Độ đen qui dẫn 1 . ị_\
Fi diện tích bề mặt bị vật bọc (vật nhỏ ) m2

F2 diện tích bề mặt vật bị bọc ( vật lớ n ) m 2
Chú ý : nếu hai tầm phẳng hoặc hai vật bọc nhau là vật trắng tuyệt đối ( vật có hệ
số phản xạ R = 1, hệ số hấp thụ A và độ đen e, A = e = 0 ) thì độ đen quy dẫn £qd = 0
và nhiệt trao đối bức xạ qi2 = 0 hay Qi2=0
2.1.2. Ví dụ
1. hai tấm phẳng đặt song song tấm thứ nhất có nhiệt độ 527°c độ đen £i =0.8 và
tấm thứ hai có nhiệt độ t2 = 27°c độ đen qui dẫn là £2 = 0.6. tính khả năng bức xạ giữa
hai tấ m :
64
do là hai tấm phẳng ta có độ đen qui dẫn là

1 1
= 0.522
. 1 1
----1------ 1 ----- 1------- 1
0.8 0.6
Lượng nhiệt ữao đổi bằng bức xạ giữa hai tẩm phẳng
T2 'ì ( 527 + 273V ( 27 + 273

?12 —£qd^ũ _ZLÌ = 0.522*5.67* = 11883*
100 100 100
w/m2)
2. Hai hình lập phương có cạnh là 5 và 200 cm. xác định độ đen quị dẫn của hệ
thống hai vật bọc nhau, tính lượng nhiệt trao đổi giữa hai vật biết độ đen qui dẫn của
hai vật lần lượt là 0.4 và 0.5 và nhiệt độ của vật 1 là 30°c và của vật hai là 15°c
Giải
Do là hai vật bọc nhau nên có thể tính độ đen theo công thức
1 1
= 0.39
1 6*0.05
- i + F' 1-1 - — -1
ex F1 0.4 6*0.22 0.5
Nhiệt lượng mà hai vật ừao đổi là:
( 30 + 273 15 +273 Ỵ
= 0.39*5.67*1.5*10-4 0.0053
Ổ l2 - £ qd c „ F \
- ì
.100 ) 100 100 100 J
Đối với cả hai bài tập này có một chương trinh chung bận có thể nhập dòng lệnh
sau đó cho chay eiải bằne trone Mathlab:
function BT2
tl=input('nhiet do cua vat 1 :');
t2=input('nhiet do cua vat 2 : ');
el= input('do den qui dan cua vat 1 :');
e2=input('do den qui dan cua vat 2: ');
colorînput('hai cat co m au:');
set=input('hai vat duoc dat nhu the nao: ');
Tl=tl+273;
T2=t2+273;
E = radiation(color,set,el,e2,Tl,T2);
disp('....................................... ');
disp('nhiet luong trao doi cua hai vat la :');
Q-E
65
Với bài 1 ta chạy chương trình và cho kết quả như sau:
» BT2
nhietdo cua vat 1: 527
nhiet do cua vat 2: 27
do den qui dan cua vat 1: 0.8
hai cat co mau: 'black'
hai vat duoc dat nhu the nao: 'songsong'
eqd =
0.5217
E =
1.1877e+004
nhỉet luong trao doi cua hai vat ỉa:

Ổ=
1.1877e+004
Vói bài tập hai ta có chạy cùng chương trình đó nhưng nhập điều kiện khác vào ta
thu được kết quả như sau:
»BT2
nhiet do cua vat 1: 30
nhiet do cua vat 2:15
hai cat co mau: 'bỉack'
hai vat duoc dat nhu the nao: 'bocnhau'
vat 1 co hình: 'cube'
nhap gia trị cua canh hmh ỉap phuong: 0.005
vat 2 co hình: 'cube'
nhap gia tri cua canh hmh lapphuong: 0.2
eqd
0.3999
nhiet luong trao doi cua hai vat la:

Ô=
66
0.0053
2.2. Truyền nhiệt và thiết bị trao đổi nhiệt
2.2.1 Tóm tắt lí thuyết

Truyền nhiệt
Truyền nhiệt qua vách phẳng
q = k(tn - t i 2) W/m2
hệ số truyền nhiệt của vách phẳng n lớp
k = ------- l— ------
i + ỷ ả + i W / m 2K
a, 1 Ắ; a 2
t n , t 12 nhiệt độ của môi chất nóng và lạnh
di ,a2 hệ số tỏa nhiệt từ bề mặt đến môi chất w/m 2 K
5; A.Í chiều dài và hệ số dẫn nhiệt của lóp thứ i
truyền nhiệt qua vách trụ
ql = kl (tu - ti2) W/m
k ,= ------------------ ỉ— -------------
1 I Ỳ 1 ln 1+1 I 1 W/mK
a xiĩdx ! lĩTẦị dị a 1ĩiả1
Trong đó kj là hệ số ừuyền nhiệt qua vách trụ n lớp
thiết bị trao đồi nhiệt
các phương trình cơ bản tính toán thiết bị trao đổi nhiệt loại vách ngăn
phương trình truyền nhiệt
Q = kFAt w
Trong đó
Q là lượng nhiệt trao đổi giữa hai môi chất, w
F diện tích bề mặt ừao đổi nhiệt, m2
k là hệ số truyền nhiệt của thiết bị ừao đổi nhiệt W/m2K
Àt độ chênh lệch nhiệt độ trung binh
phương trình cân bằng nhiệt
Q= - 1n = G2Cp2(t2’ - 12”) w
Chỉ số 1 la của chất lỏng nóng còn chỉ số 2 là của chất lỏng lạnh
Các kí hiệu “ ‘ “ các thông số đi vào thiết bị
Các kí hiệu “ “ “ các thông số đi ra thiết bị
G lưu lượng kg/s
67
G = vp
V là lưu lượng thể tích m3/s
p là khối lương riêng của chất lỏng kg/m3
Cp nhiệt dung riêng khối lượng đẳng áp J/kgK
Q = Wi 5ti = w 2 5t2 ,w , = G ,cpl; w 2 = G2Cp2
độ chênh lệch nhiệt độ trung bình ĩogarit của thiết bị cùng chiều và ngược
chiều
1n-=í-
At1
Đối với chất lỏng song song cùng chiều
Ati —t]’ —\à
At2 = t]” - t 2”
Đối với các chất lỏng song song ngược chiều
Àtj = t]’ —Í2”
At2 = t r - V
2.2. xác định diện tích bề mặt thiết bị trao đổi nhiệt
F = kAt
£
2.1.2. Ví dụ
3. Vách phẳng hai lóp chúng đều có chiều dày là 0.25m hế số dẫn nhiệt của hai
lớp lần lượt là 0.348 W/mK và 0.695 W/mK hệ số tỏa nhiệt tương ứng của chủng là
34.8 W/m2K và 11.6W/m2K nhiệt độ lớp trong cùng là 1300°c nhiệt độ của lóp ngoài
cùng là 30°c. tính mật độ dòng nhiệt
Lời 2Ìải
Ta có do vách là vách phang nên ta dùng côn thức sau tính hệ số tỏa nhiệt của hệ
ỵ- _______ Ị__________________Ị__________ _ 0 838
I +ỷ i +l 1 | 0.25 | 0.25 |1' W/m2K
ax Ỷ Ậ «2 34-8 0 348 0 695n -6
Mật độ dòng nhiệt
Q = kAừ=0.838(1300 - 30)=1064 w /m 2
4. một ống dẫn hơi làm bằng thép đường kính là 200/216 mm. hệ số dẫn nhiệt
bằng 46W/mK. Được bọc bằng một lớp bảo cách nhiệt dày 120 mm có hệ số dẫn
nhiệt là 0.116W/mK. Nhiệt độ hơi bằng 300°c. hệ số tỏa nhiệt từ bề mặt ừong tới
không khí là 116W/m2K và hệ số tỏa nhiệt tù bề mặt ngoài tới không khí xung quanh
là 10 W/m2K. xác định mật độ dòng nhiệt
Giải
68
Do là vách trụ nên áp dụng phương trình sau:

d 3 = d2 +25=216+2.0.12=0.456(m)
k t= _______ _ _ i ___________ = _____________________ 1____________________
' —-— 1 + V- 1 ,to
> —— _ 4^ ±L
+i + — -—
1 ---- 1---- + — -—
1 ln:^
1 216
—+ ----- 1---- ln——
, 456 + ----- 1-----
a xJtdx di a 1nd1 116^0.2 2^216 200 2^0.116 216 10^0.456
(W/mK)
Q=kÁt=0.9*(300-25)=247.5(W/m)
Với cả hai bài này ta có một chương trình chung nhập dòng lệnh sau đó cho chạy
trên mathlab ta sẽ có kết quả của cả hai bài tập
íiinction BT3
tl=input('nhiet do cua lop trong cung:');
t2=input('nhiet do cua lop ngoai cung:');
lamda=input('nhap duoi dang ma ừan he so dan nhiet cua cac lo p :');
al=input('nhap he so toa nhiet cua lop trong cung:');
a2=mput('nhap he so tao nhiet cua lop ngoai cung:');
delta=input('nhap duoi dang ma tran gia tri cua chieu day cac lo p :');
vach=input('loai vach ma nhiet truyen q u a:');
k=heat(vach,tl ,t2,al ,a2,lamda,delta)
dispC— — .......
disp('nhiet trayen qua vach la :')
Q=k*(tl -t2)
Với BT3 cho chạy thu được kết quả sau :
BT3
»BT3
nhiet do cua ỉop trong cung: 1300
nhiet do cua lop ngoai cung: 30
nhap duoi dang ma tran he so dan nhiet cua cac ỉop: [0.348 0.695]
nhap he so toa nhiet cua lop trong cung: 34.8
nhap he so tao nhiet cua lop ngoai cung: 11.6
nhap duoi dang ma tran gia tri cua chỉeu day cac lop: [0.25 0.25]
loai vach ma nhiet truyen qua: 'phcmg'
q=
1.0645e+003
Với BT4 có kết quả như sau
69
»BT3
nhiet do cua lop trong cung: 300
nhiet do cua lop ngoai cung: 25
nhap duoi dang ma tran he so dan nhiet cua cac lop: [216 0.116]
nhap he so toa nhiet cua lop trong cung: 116
nhap he so tao nhiet cua lop ngoai cung: 10
nhap duoi dang ma tran gia tri cua chieu day cac lop: [0.008 0.12]
loai vach ma nhiet truyen qua: 'ừĩi'
nhap gia tri cua duong kinh ong trong cung: 0.2
k=
0.9019
nhiet truyen qua vach la:

Q=
248.0206
5. Trong một thiết bị trao đổi nhiệt cần làm nguội 275Kg/h chất lỏng từ 120°c tới
50°c chất lỏng nóng có nhiệt dung riêng 3.04 kJ/kgK chất lỏng lạnh có lưu lượng
lOOOkg/h nhiệt độ vào thiết bị là 10 nhiệt dung riêng
c2“ 4.18 kJ/KgK. Biết hệ số truyền nhiệt thiết bị là k=1160W/m2K. tính diện tích
của truyền nhiệt của thiết bị trong các trường hợp :
chất lỏng chuyển động song song cùng chiểu
chất lỏng chuyên động song song ngược chiều
Lời giải
nhiệt lượng do chất lỏng nóng thải ra
Q= GiCiÀt’= 275/3600.3.04.10"3.(120 - 50) = 16255.5 w
Nhiệt độ của chất lỏng lạnh ra được xác định từ pt cân bằng nhiệt
_ 10| 275*3.04*(120-50)_ 21„c
1000*4.18
Độ chênh lệch nhiệt độ
c)song song cùng chiều
Ati = t i’ - t 2’ = 1 2 0 - 10 = 110
Át2 = t j - 12” = 5 0 - 2 4 =26
Át, 26
70
= 0.24m2
Q _ 16255
Diện tích F= 0.22m 2
kAt 1160*64
function BT4
Gl=input('nhap luu luong cua dong nong:');
G2=input('nhap luu luong cua dong lanh:');
tnv=input('nhap nhiet do ban dau cua dong nong:');
tnr=input('nhap nhiet do di ra cua dong nong: ');
tlv=input('nhap nhiet do ban dau cua dong lanh:');
tlr=input('nhap nhiet do dong lanh di r a :');
cl=input('nhap nhiet dung rieng dang ap cua dong nong cua: ');
c2=mput('nhap nhiet dung rieng dang ap cua dong lanh: ');
k=input('he so truyen nhiet cua chat long:');
moving=input('hai dong chuyen dong ra sao:');
deltatn=(tnv-tnr);
Q=heatquatity(cl ,G1 ,deltatn);
x=heatbalance(Gl ,G2,cl,c2,tnv,tnr,tlv,tlr);
t=average(moving,tnv,tnr,tlv,x);
disp('
disp('dien tich be mat truyen nhiet la :');
F=surface(Q,k,t)
Kết quả thu được là:
»BT4
nhap ỉuu luong cua dong nong: 275/3600
nhap ỉuu luong cua dong lanh: 1000/3600
nhap nhiet do ban dau cua dong nong: 120
71
nhap nhiet do di ra cua dong nong: 50

nhap nhiet do ban dau cua dong lanh: 10
nhap nhiet do dong lanh di ra: 'no ’
nhap nhiet dung rieng dang ap cua dong nong cua: 3.04e3
nhap nhiet dung rieng dang ap cua dong lanh: 4.18e3
he so truyen nhiet cua chat long: 1160
hai dong chuyen dong ra sao: 'comoving'
Q=
1.6256e+004
nhiet do cua lop chat long lanh di ra la:
x =
24
t=
58.2369
dien tich be mat truyen nhiet la:

F=
0.2406
Với phần b kết quả là chú ý khi nhập thay đổi hai dòng chuyển dộng ra sao nhập
‘reverse’(ngược chiểu)
»BT4
nhap luu luong cua dong nong: 275/3600
nhap luu luong cua dong lanh: 1000/3600
nhap nhiet do ban dau cua dong nong: 120
nhap nhiet do di ra cua dong nong: 50
nhap nhiet do ban dau cua dong lanh: 10
nhap nhiet do dong lanh di ra: 'no'
nhap nhiet dung rieng dang ap cua dong nong cua: 3.04e3
nhap nhiet dung rieng dang ap cua dong lanh: 4.18e3
he so truyen nhiet cua chat long: 1160
hai dong chuyen dong ra sao: 'reverse'
Q=
1.6256e+004
72
nhỉet do cua ỉop chat long lanh di ra la:

X -
24
t=
63.9657
dien tich be mat truyen nhiet la:

F=
0.2191
73
Chương 3: Kỹ thuật tách chất
Tách chất là một trong những giai đoạn rất quan trọng ừong quy trình công nghệ
hóa học. Dựa vào các đặc tính hóa lí khác nhau của các chất mà người ta tách riêng các
chất ra khỏi hỗn hợp. Trong thực tế có rất nhiều phương pháp được sử dụng, một trong
số đó là phương pháp chưng cất dựa vào sự phân biệt nhau về nhiệt độ sôi của các cấu
tử ừong hỗn hợp lỏng cần tách.
3.1 Cân bằng lỏng hoi
3.1.1 Tóm tắt ỉỷ thuyết

-) Định luật Raoult’s:
Pi = x rPũi (3.1)
-) Định luật Dalton
p, =y,-P (3.2)
-) Ket hợp hai định luật trên cho hệ hai cấu tử ta có:
p-pn
1”
p 01 Pai
(3.3)
(3.4)
(3.5)
(3.6)
H ay :
(3.7)
F=v +L với từng thành phần thì: ZịF = XịL + y ỵ (3.8)

Hay z i =Ợ C i - \ ) x i v '+ x ỉ ( 3 .9 )
74
(3.10)
(^ -iy + 1
Trong đó :
p , : áp suất hơi riêng phần của cấu tử i
X. : phần mol của cấu t ử I trong hỗn hợp
p m: áp suất hơi bão hòa của cấu tử i ở trạng thái
nguyên chất
P: áp suất hơi bão hòa của hệ
A,B,C : các hằng số
ớ: nhiệt độ (°C)
z,: Phần mol của nguyên liệu vào
F: Tổng lưu lượng mol
V: Lưu lượng mol của pha hơi
L: Lưu lượng mol của pha lỏng
3.1.2 Ví dụ
1. Cho các hằng số của phương trình Antoine
Với benzen: kl = 6.90565 k2 = 1211.033 k3 = 220.79
Với toluen: kl = 6.95334 k2 = 1343.943 k3 = 219.377
Trong đó p a(mm Hg)
Xác định thành phần pha hơi trong hỗn họp (benzene và toluene) biết rằng thành
phần của chúng trong pha lỏng là 0.5 benzen và 0.5 toluen ở 65 ° c . Hỗn hợp được bay
hơi ở áp suất 1atm.
Lời siải
Áp dụng phương trình (3.4)
Cho Benzen:
IgtPoi) = 6.90565 - [1211.033/(65 + 220.70)] = 2.668157
<=> p n = 465.75 mm Hg hay 62.10 kN/m2
Cho Toluen :
lgCPo2) = 6.95334-[1343.943/(65 + 219.377)1 = 2.22742
<=> p 02= 168.82 mtn Hghay 22.5 kN/m2
=> ^ = (0 .5 0 X62.10) = 31.05kN/m2
75
p 2= (0.50 X22.51) = 11.255 kN/m2

p = p l +p1 = 42.035 kN/m2
Sử dụng phương trình (3.2):
yl - (31.05/42.305) = 0.734
Và:
y 2 = (11.255/42.305) = 0.266
2. Cho một hỗn họp lỏng ở 46.25 ° c và 4.8253 Bar có thành phần dòng vào như
trong bảng 3.1. Tính thành phần mỗi chất ừong pha lỏng (x.) và thành phần trong pha
h ơi (>»,).
Bảng 3.1:Thảnh phần phần mol và giá ừị cân bằng lỏng hơi của các cấu tử
Cấu tử Phần mol(z ) Giá trị cân bằng pha(iỉrí)
Propan 0.1 6.8
n-butan 0.3 2.2
n-pentan 0.4 0.8
n-octan 0.2 0.052
lời 2Ỉải
Kết hợp phương ừình (3.7) và phương trình (3.10) ta có:
y =0
Thay số:
( 6 .8 - l) x 0 .1 ( 2 .2 - l) x 0 .3 ( 0 .8 - l) x 0 .4 (0.0 5 2 -1 ) X 0.2
l + (6 .8 -l)v ' l + (2 .2 -l)v ' l + (0 .8 -l)v ' 1 + (0 .0 5 2 -l) v '”
<=> 0.58 0.36 0.08 0.1896

l + 5.8v' 1 + 1.2v' l-0 .2 v ' 1+ 0.948v'_
<=> v'= 0.4258
Thay số vào (5.5) và (5.10) ta có bảng giá trị
76
Bảns 3.2 Kết quả

Cấu tít Thành phần pha lỏng( X ị ) Thành phẩn pha hơi( y, )
Propan 0.028 0.196
n-butan 0.1985 0.4368

n-pentan 0.4372 0.3498
n-octan 0.3354 0.0174
Giải bằne Matlab

Khởi động Matlab rồi gõ » e d i t
Màn hình sẽ hiện lên cửa sổ soạn thảo
Với ví dụ 1 chúng ta có thể dùng code như sau:
fonction btl
kl = input('cac hang trong phuong ừinh Antoine cho chat thu nhat:');
%phải nhập dưới dạng ma trận
k2 = input('cac hang toong phuong trinh Antoine cho chat thu hai:');
%phải nhập dưới dạng ma trận
t = inputCNhiet do cua he(do C):');
%nhiệt độ phải đổi ra độ c
p l = 10Aantoine(kl,t);
p2 = 10Aantoine(k2,t);
p = p l +p2;
yl = pl/p;
y2 = p2/p;
disp(—--......................................................... ')
fprintf('Thanh phan cua chat thu nhat trong pha hoi la:yl =
%f(%%)\n',yl * 100)
fprintfi'Thanh phan cua chat thu hai trong pha hoi la:y2 =
%f(%%)\n',y2* 100)
chủng ta lưu file với tên ‘ btl .m ’
tuy nhiên muốn cho chương trình có thể chạy được thì chúng ta phải có hàm
‘antoine.m’.
Chú ý: Cả hàm ‘antoine.m’ lẫn file ‘btl .m’ đều phải nằm trong Current directory,
code:
function y = antoine(k,t)
77
y = k(l)-(k(2)/(k(3)+t));
return
Để giải ví dụ 1, tù cửa sổ CommandWindow của Matlab gõ lệnh
» b tl
cac hang trong phuong ừmh Antoine cho chat thu nhat:[6.90565,1211.033
,220.79]
cac hang trong phuong trình Antoine cho chat thu hai:[6.95334, 1343.943,
219.377]
Nhiet do cua he(do C):65
Thanh phan cua chat thu nhat trong pha hoi la:yl = 73.396546(%)
Thanh phan cua chat thu hai trong pha hoi la:y2 = 26.603454(%)
Vói ví dụ 2 chúng ta có thể gõ code như sau:
function bt2
global z K
z = input('phan mol cua nguyen vao (zi) la:');
%cho thành phần mol của các cấu tử trong hỗn họp vào dưới dạng
%ma trận
K = input('Gia tri can bang pha (Ki) là:');
%cho giá trị cân bằng pha của các cấu tử đó dưới dạng ma ừận
disp('..................................................... ')
V = fzero(@vpequil,0.2);
%hàm vpequil(v) có thể tim hiểu ừong 7.5.9
n = length(z);
for i = 1:n
x(i) = z(i)/((K(i)-l)*v+l);
y(i) = K(i)*x(i);
end
%cho kết quả phần mol của các cấu tử ữong pha lỏng
disp('thanh phan mol cac chat trong pha long la:')
X
%cho kết quả phần mol của các cấu tử ữong pha hơi
disp('thanh phan mol cac chat trong pha hoi la:')
y
chúng ta lưu file với tên ‘ bt2.m ’
78
cững giống như ví dụ 1, muốn cho chương trình có thể chạy được thì chúng ta
phải có hàm ‘vpequil.m’
code:
function y=vpequil(v)
global z K
sum-0;
n = length(K);
for i=l:n
tu = (K(i)-l)*z(i);
mau= l+(K(i)-l)*v;
sum = sum + tu/mau;
end
y=sum;
Để giải bài tập 2 , tù cửa sổ CommandWindow của Matlab gõ lệnh
» bt2
phan moỉ cua nguyen vao (zi) la:[0.1 0.3 0.4 0.2]
Gia ừi can bang pha (Kỉ) ỉà:[6.8 2.2 0.8 0.052]
thanh phan mol cac chat trong pha long la:

X
0.0288 0.1985 0.4372 0.3354

thanh phan moỉ cac chat trong pha hoi ỉa:
y=
0.1960 0.4368 0.3498 0.0174
Ket quả hoàn toàn giống như trong bảng 5.2
3.1.3. Bài tập

1 .Một dòng hỗn họp chảy vào thiết bị ở 45.2650c và 6.8932 Bar. Xác định thành
phần mỗi chất trong pha lỏng ( x j và thành phần trong pha hơi (}>,.).
79
Bànz 3.3 Thành phần phần mol và giá ừị cân bằng lỏng hơi của các cấu tử
Chất lbmol/h K,
Propan 20 3.7
n-butan 30 1.4
n-pentan 50 0.6
llbmol = 453.6 mol
2. Cho các hằng số của phương trình Antoine

Với benzen: kl = 6.90565 k2 = 1211.033 k3 = 220.79
Với toluen: kl = 6.95334 k2 = 1343.943 k3 = 219.377
Trong đó p 0(mm Hg)
Xác định thành phần pha hơi và pha lỏng ừong hỗn hợp (benzene và toluene) biết
rằng áp suất chung của cả hệ là 56 kN/m2 ở 950c .
3.2 Phương pháp McCabe Thiele
3.2.1 Tóm tắt lý thuyết

-) Phương trình cân bằng vật chất
F=D+B (3.11)
F .x f = D .x d + B .xb ( 3 .1 2 )
Trong đó :
F: lưu lượng nguyên liệu đưa vào cột cất (mol/h)
D: lượng trong phần cất (mol/h)
B: lượng chất trong phần đáy cột cất
XF : Phần moi của thành phần A trong nguyên liệu vào
XD : Phần mol của thành phần A trong phần cất
XB : Phần moi của thành phần A trong phàn luyện ở đáy
-) Phần trên cột cất(phương trinh đường chưng cất)
y ^ ỉ = ệ xM + yxD (3-13)
Trong đó:
L: lưu lượng mol lỏng trong toàn dòng từng giai đoạn (mol/h)
V: lưu lượng mol hơi trong toàn dòng từng giai đoạn (mol/h)
80
D: Sản phẩm phần cất (mol/h)

XD: Phần mol của thành phần dễ bay hơi hơn trong phần cất
xn:Phần mol của thành phần dễ bay hơi hơn ừong phần lỏng lấy ra ở giai
đoạn n
y n+Ị: Phần mol của thành phần dễ bay hơi hơn trong phần hơi lấy ra ở
giai đoạn n
Đặt R = — là tỷ số hồi lưu, phưong trình (3.13) có thể viết lại:
R 1 /'t 1
V . = *„ + _ — x n (3 .1 4 )
"+1 R + 1 ” R +1 D v ’
-) Phần dưới cột cất (phương trình đường luyện)
ym+i = ỹ xn > - ỹxB (3-15)
Trong đó:
L : toàn bộ lượng pha lỏng trong phần dưới cột cất (mol/h)
V : dòng hơi trong từng giai đoạn ở phần dưới cột cất
Do V =L —B nên (5.15) có thể viết lạ i:
y”*1=TZ~ãXm~TZ~ãXB
L —D L —Đ (3-16)
Sử dụng cân bằng vật chất ta có:
F + L +V =V + L (3.17)
L-L
Đăt q =------ = (nhiêt đê chuyên hóa 1 mol nguyên liêu thành hơi bão hòa)/(nhiêt
F
hóa hơi 1 mol nguyên liệu), phương trinh (3.16) có thể viết:
= Lr +
ỉ q+FF -- B
B X»~ L-\-qF-B
r-L. (3-18)
Phương ừĩnh đường q ửiu được:
y =— x- —T
q- 1 q- 1
(3-19)
Giá trị của q:
Điều kiên của nguyên liêu Tính giá ừi a
Lỏng quá lạnh q = l + CPL(tB- t F)/Ầ
Tại điểm sôi q= 1
Một phần hơi q =f L
Điểm ngưng tụ q=0
81
Trong đó:
CPL :nhiệt dung dạng lỏng
cpv:nhiệt dung dạng hơi
tp : nhiệt độ của nguyên liệu vào
tD: điểm ngưng tụ
tB : điểm sôi
f L : phần mol dạng lỏng
Ằ :nhiệt hóa hơi cho 1 mol
-) Lượng hồi lưu tối thiểu

xD ~ y P
m —- (3.20)
Hay r^ = 7Z~ (3.21)

l-m
-) Sự hồi lưu toàn phần
Tính số đĩa lý thuyết
X D (1 — X B )
lg
XB ặ - X D )
N*>=- (3.22)
!g«av
Trong đó :
a - y( ỉ - x ) (3.23)
x ( y - 1)
«av 'ị& to p & bottom (3.24)
-) Cân bằng năng lượng ttong cột cất

F.hF +Qb = D.hD + B.hB + Qc (3.25)
Trong đó: Qc: công suất bộ ngưng tụ
Qb : công suất nồi hơi
h: các entapy
Nếu ngưng tụ tại điểm sôi của thì:
Qc =DịR + \)1 (3.26)
Trong đó: Ả : nhiệt ngưng tụ
82
Nêu tf ờ điều kiện chuẩn, nguyên liệu 100% là lỏng, và tD,tBlà nhiệt độ của sản
phẩm trên cột cất và đáy cột cất
hF =0, hD = CpD.(tD— và hs = CpB(¿2 ~tp)
Trong đó : CPDi c ps: nhiệt dung riêng đẳng áp của phần trên và phần đáy
cột cất
Qb = D.Cpd(tD—tỵ) + B.Cpß (tB—tỵ)+Qc (3 .27)
3.2.2 Vỉ dụ
3.Một hỗn hợp chứa 50% mol n-pentan yà 50% mol n-hexan,sản phẩm thu được
ở phần trên cột cất chứa 95% mol pentan và sản phẩm ớ đáy cột chứa 5% mol pentan.
Chưng cất ừong ớ áp suất latm. Giản đồ cân bằng như trong hình 3.1. Nguyên liệu
đưa vào đang ở điểm sôi.
Hình 3.1
ã) Tính tỷ số hồi lưu nhỏ nhất
b) Tính sổ đĩa lý thuyết nhỏ nhất
4.Một cột cất phân đoạn được thiết kế để tách 500 Kg/h hỗn hợp gầm 40% mol
benzene và 60% mol toluene, sản phẩm trong phần cất gồm 97% mol benzene và sản
phẩm ở đáy chứa 98% mol toluene. Quá trình chưng cất được thục hiện ở áp suất khí
quyển. Cho biết:
Nhiệt độ sôi của nguyên liệu lỏng = 368 K
Nhiệt độ của nguyên liệu đầu vào = 293 K, nhiệt dung riêng nguyên liệu là :
1.6xl02 KJ/(kmoLK)
Nhiệt hóa hơi của benzene = 3.08X104kJ/kmol
Nhiệt hóa hơi của toluene = 3.33 X104kJ/kmol
Giản đồ cân bằng của hệ benzene- toluen như trong hình 3.2
ọ
I
0
1
is
0.2 0.4 0.6 08
X, mol ừaction benzene in lỉquid

0 9 /0 9 /0 0 16:20
Hình 3.2
a) Tính hệ số góc của đường q
b) Tỷ số hồi lưu thực tế lớn gấp 1.5 lần . Hãy xây dựng cân bằng vật chất,
công suất bộ chưng cất,công suất tháp chưng cất Nhiệt độ phần trên cột cất = 355K,
nhiệt độ phần đáy cột cất = 382.8K . Bỏ qua sự thất thoát nhiệt.
c) Tính lưu lượng lỏng và hơi trên toàn dòng ở phần ừên và đáy cột cất
Lởi ẹiải
3. a)Dựa vào hình 3.1 ta thấy:
—^ = 0.542
iL +1
<=> 2L = 0.75
b)Từ hình 5,1 ta thấy: XD = 0,95 , thành phần cân bằng hơi yD = 0.98
XB =0.05 , thành phần cân bàng hơi y B = 0.11
Áp dụng phương trình (3.23) ta có:
0.98 1-0.95
= 2.58
1-0.98 0.95
0.11 1-0.05
a bottom - x--- = 2.35
1-0.11 0.05
a av = V2.58x2.35 = 2.46
Theo phương trình Fenske:
84
/ 0.95x(l-0.05)^
log
=— V0.05
--- XV(1-0.95)
- -1=5.5
log(2.46)
4. a) Hỗn hợp nguyên liệu vào có nhiệt độ 293 K nên có thể coi là điều kiện
nguyên liệu lỏng lạnh nhiệt. Do đó để tăng 1 mol hỗn họp vào lên 368K =
1X(1.6 X10* )(368 -293)=12000 kJ/kg.mol
Ẳ của hỗn hợp = 0.4(3.08xl04) + 0.6(3.33xl04) =3.23xl04 kJ/kmol
1 2 0 0 0 +3.23xio4 , „„
--------------------- = 1.37
3.23x10
Hệ số góc của đường q = —— = = 3.7

6 q - 1 1.37-1
b) Vẽ đường q với hệ số góc là 3.7 đi qua điểm ( xF,yF) cắt đường cong tại p.
Nối (xD,yD) vói p. Ta thấy đường này cắt trục tung tại (0, 0.425)
Do đó : = 0.425 <=>Rm=1.28
Rm+1
m
Tỷ số hồi lưu thực là: 1.5x1.28 = 1.92
Khối lượng moi của nguyên liệu = 0.4(78) + 0.6(92) = 86.4
Lưu lượng mol nguyên liệu = 5000/86.4 = 57.87 kmol/h , XD = 0.97 , XB = 0.02
Thay số vào (3.11) và (3.12)
[D+B =57.87 „
_ _ <=>D = 23.16 B = 34.71 kmol/h
(0.97Z) + 0.025 = 0.4 X57.87
Công suất bộ ngưng tụ:
Ẳ = 0.97(3.08 xio4) + 0.03(3.33 xio4) = 3.09xl04 kJ/kmol
Qc =23.16(1.92 + 1)(3.09x104) = 2.09x104 kj/h
Công suất tháp chưng cất:
tD=355K , tB= 382.8.sr
Thay số vào phương ừình (3.27)
Qb = 23.16 X1.6 X102(355 - 293)+ 34.71x1.6 xio2(382.8-293) +2.09 xio6 = 2.82xl06
kj/h
c) L = RD = 1.92(23.16) = 44.47 kmol/h
V = D(R+1) = 23.16(1.92+1) = 67.63 kmol/h
L = L +qF = 44.47 + 1.37(57.87) = 123.75 kmol/h
V = V - (1- q)F = 67.63 - (1 - 1.37)(57.87) = 89.04 kmol/h
Giải bằns Matlab
85
3 .Với ví dụ 3 chúng ta có thể viết mã code :

function bt3
x f = input('phan moi cua chat trong nguyen lieu ( xf
xd = input('phan mol cua chat trong phan ừen cot cat ( xd
yd = input('phan mol cua thanh phan hoi ( yd
xb = input('phan mol cua chat trong phan day cot cat ( xb
yb = input('phan mol cua thanh phan hoi ( yb
yO = input('tung do cua diem cat teen gian do:');
%tính R min
Rm = xd/yO-1;
%tính số đĩa lý thuyết theo phương trình Fenske
Nm = Fenske(xd,yd,xb,yb);
disp('..................... kq..........................')
fprintf('Ty so hoi luu toi thieu la: %f\n',Rm)
íprintf('So di ly thuyet tinh theo phuong trinh Fenske la: %f\n',Nm)
lưu file với tên : ‘bt3.m’ và chú ý như các ví dụ trước. Muốn file chạy được thì
phải có hàm ‘Fenske.m’(sẽ được tìm hiểu hàm này trong chương 7) tồn tại trong
CurrentDirectory. Để giải ví dụ 3,ừên của sổ CommandWindow gõ lệnh » b t3 và
nhập số liệu theo chỉ dẫn ta được :
» bt3
phan mol cua chat trong nguyen lieu ( x f ):0.5
phan mol cua chat trong phan tren cot cat ( xá ):0.95
phan mol cua thanh phan hoi ( yd):0.98
phan mol cua chat trong phan day cot cat ( xb ):0.05
phan mol cua thanh phan hoi ( yb ):0.11
tung do cua diem cat tren gian do:0.542
---------------- kq-------------------
Ty so hoi ỉuu toi thieu ỉa: 0.752768
So di ly thuyet tinh theophuong trinh Fenske ỉa: 5.539277
4. Còn với ví dụ 4
Code:
function bt4
global Xhv Cp R
86
%lamda:nhiệt hóa hơi của 1 mol hỗn hợp

%Cp: nhiệt dung của hỗn hợp, có thể là Cpl hoặc Cpv tùy điều kiện của
nguyên liệu đầu
%Cpl:nhiệt dung của hỗn hợp lỏng
%Cpv:nhiệt dung của hỗn hợp ở pha hơi
X = input('nhap ma tran phan moi cua cac thanh phan trong hon hop:');
hv = input('nhap ma tran nhiet hoa hoi tuong ung cua cac thanh phan:');
Cp = input('nhiet dung cua hon hop:');
t = input('tB(trong truong hop long qua lanh) hoac tD(ừong truong hop
nhiet qua soi):');
tf = input('nhiet do cua nguyen lieu vao(tí):');
yO = input(toa do cua duong chung cat voi truc tung:');
F = input('luu luong tong cong cua nguyen lieu vao(kg/h):');
xd - mput('phan mol cua benzen trong cot cat(xd):');
td = inputCnhiet do cua cot cat:');
xb = input('phan mol cua benzen trong phan day cot cat(xb):');
tb - inputCnhiet do cua phan day cot cat:');
dk = 'long qua lanh';
%tính R min
Rm = xd/yO-1;
%tính R thực tế
R = 1.5*Rm;
%tmh hệ số góc của đường q
[q, q_slope] = q_line(dk,t,tf);
m = x(l)*78 + x(2)*92;
Fn = F/m;
%tính D và B
[D,B] = materìal(Fn,x(l),xd,xb);
X = [xd 1-xd];
%tính Qc và Qb
[Qc,Qb] = Heat(D,td,B,tb,tf);
L = R*D;
V = D*(R+1);
Ltb = L + q*Fn;
vtb = V - (1 - q)*Fn;
disp('.......................kq...........................')
fprintf('he so goc cua duong q la: %f\n',q_slope)
fprintf('D = %f (kmol/h)\n',D)
fprintf('B = %f (kmol/h)\n',B)
fprintf('cong suat cua bo ngung tu la: %f (kJ/h)\n',Qc)
fprintf('cong suat cua thap chung cat la: % f (kJ/h)\n',Qb)
fprintf('L = %f (kmol/h)\n',L)
fprintf('V = %f (kmol/h)\n',V)
fprintf('Vtb = %f (kmol/h)\n',Vtb)
fprintf('Ltb = %f (kmol/h)\n’,Ltb)
Lưu với tên ‘bt4.m \ Giải ví dụ 4 ta gõ vào CommandWindow
» bt4
nhap ma franphan mol cua cac thanh phan trong hon hop:[0.4 0.6]
nhap ma tran nhiet hoa hoi tuong ung cua cac thanh phan:[3.08*10*4
3.33*10A4]
nhietdung cua hon hop:1.6*10A2
tBậrong ừuong hop long qua lanh) hoac tD(trong ừuong hop nhiet qua soi):368
nhiet do cua nguyen lieu vao(tj):293
toa do cua ãuong chung cat voi truc tung:0.425
ỉutt luong tong cong cua nguyen lỉeu vao(kg/h):5000
phan mol cua bemen trong cot cat(xd):0.97
nhietdo cua cot cat:355
phan mol cua bernen trong phan day cot cat(xb):0.02
nhiet do cua phan day cot cat:382.8
----------------- kq--------------------
he so goc cua duong q la: 3.691667
D = 23.148148 (kmol/h)
B = 34.722222 (kmol/h)
cong suat cua bo ngung tu la: 2089443.763617 (kJ/h)
cong suat cua thap chung cat la: 2817962.282135 (kJ/h)
L = 44.526144 (kmol/h)
v= 67.674292 (kmol/h)
Vtb = 89.174120 (kmoỉ/h)
Ltb = 123.896342 (kmol/h)
88
Chương 4: K ỹ thuật phản ứng

4.1 Nhiệt động học
4.1.1 Tóm tẳt lý thuyết

A H ^ỵA H ^-ỵA H ^ (4.1)
AHr = AH Ta + £ ACPdT (4.2)
Trong đó:
ACp = ni ,Cpi (product) - ^ nt ,CPi ('react)
AG = -R T h íK (4.3)
Định luật Van’t H o ff:
d \n K ah"
(4.4)
dT ~ RT 2
Neu AHị thay đổi không đáng kể theo nhiệt độ thì định luật trên còn được viết
1 1
dưới dạng: ln— = (4.5)
Kt R T
V1! T
Ẳu
Tốc độ phản ứng :
(4.6)
V dt
Bảns 4.1 : Tốc độ của các phản ứng đơn giản
phản ứng bậc tốc độ dạng tích phân t(bánhuỷ)
A -> p 0 = CM - C À =kt *1/2= —

dt 2k
A -> p 1 z^ L = kCÁ ln^ - = kt í1/2 = —ln2

dt A CA k
2A -> p 2 Z^£ l = ìccI -ị — - ? - = ** ty , = - ^ —

dt CA CA0 kCM
-) Phản ứng thuận nghịch

A <=> B
= ^ J - = k c A- k1CB (4.7)
dt dt 1 Á 1 B v ’
Dạng tích phân :
Cá Ca° = e k*‘ (4.8)

cÂO - C^Ẩồ
89
Trong đó :
kị (K +1) f _
c BO T+ c AO k.
kg —- K = ^~
K K +\ k,
-) Phản ứng song song
A ^B
dc. _ dCR , dCc
= (kl +k1)CẢ ; - ^ f - = kíCẨ = k1CÀ
dt dt dt
với điều kiện tại t = 0, CẨ = CẦữ;CB =Cc = 0,
K 1 -^ K - ( l- x < ) ; x c = - ^ - =
CẨữ kị + k2 ^ CM) kị + k. CA0 kị + k1 kị + k.
-) Phản ứng nối tiếp

kị k2
Â -ỳB -ỳC
Điều kiện ban đầu : tại t = 0, CA = CÃ0, CB = Cc = 0
dCẮ =]Í Q ■ đCB _ k c _ ỵ c . dCc —k2CB
dt lẨ dt 1 À 2 B dt
viết dưới dạng tích phân :
e-h‘ e-h‘ s
c.
'■'Ã =c. k *>' ’■ c. ~
—'âoíT =c. k Cc —CẨ0
1^2 ẨúỊ
Nồng độ của B lớn nhất tại thời điểm :
ln(fr2 /¿ị)
t= vàX ——
CE k l
k2 —kị \ KU
Phương ừình tốc độ biểu diễn qua độ chuyển hoá :
-) Bậc n h ất: A -> p
dX,
= k ( \ - X Ả) (4.9)
dt
Hay - \ n ặ - X Ả) = kt t ạ i t = 0 , X Ã =Q (4.10)
-) Phản ứng bậc 2 : A + B -> p
CM d^ - = kC1Aa{ \ - X i ) { M - X A) M =^
CM
C4.ll)
90
+) Neu : M *■1 thì ln ——— = CM{M - 1)kt (4.12)

’ M( \ - X A) AữK ’
+) Nếu : M = 1 thì = kCẦ0(1 —X A)2 (4.13)

dt
Hay —^----- — — = kt ,tại t = 0 thì X , = 0 (4.14)

J CM l - X A
+) Tương tụ với bậc n :
dx
- ^ = kC»Ẩ-l( l - X Ẩy (4.15)
at
Tại t = 0 có X Â = 0 thì (6.15) có thể viết dưới dạng:
(l - X Af n - \ = { n - \) C ^ k t (4.16)
k2
-) Phản ứng song song: A —} B Ẩ-ỳ,C
1n(l- X Ầ) = {kl +k1)t (4.17)
-) Phản ứng thuận nghịch:
dXA _ kị (M +1)
(4.18)
dt ~ M +x
Viết dưới dạng tích phân:
M +ì
-ln 1 - ^ kịt (4.19)
X M +X t
Trong đó X Ae là độ chuyển hoá lúc cân bằng
Với phản ứng thuận nghịch bậc 2: A + B <=> R + s
d^ - = klc 1AŨ{ \ - x ẢÝ -
dt kị
Dưới dạng tích phân:
In X ụ -V *A e -1 )X ạ —2klCẨ
X Àe - X À A y
Phương trình Michaelis - Menton :
dCA CA
- r t = - —— = kl (4.20)
dt 1+ k1CA
4.1.2 Ví dụ
1. Etyl benzene phân huỷ theo phản ứng:
C,H5C1H ỉ ^ C 6H 5C1H ì + H 1
với hằng số tốc độ phản ứng ở:
91
nhiệt độ °c ẮrxlO4
540 1.6
550 2.8
Tính năng lượng hoạt hoá của phản ứng.
2. cho phản ứng nối tiếp
A -> B -> p coi phảnứng là phản ứng bậc 1 vớikị = 1h~l và k2 =2h~l , tại thời
điểm t = 0 , Ca = 5 mol/1, Cb = Cc = 0 mol/1.Hãy biểu diễn sự thay đổi nồng độ mỗi
chất vào thời gian.
lời ziải
1. Từ phương trình Arrheneous : k = ka.e E/RT
k e~EIRT>
k, e~EIST'
lấy logarit 2 vế ta có:
ln— =
kị RT2 + R ĩ ;
JỊT Ị_ 2.8x10
<=>E = ln — = ln 8.314X 823x813 1= 311308 (J/mol)
kx Ti-Ti v1.6 x l 0 J 823-813 1

1. Tạo màn hình soạn thảo tương tự những ví dụ của chương trước. Sau đó gõ mã
code :
íunctìon btl
%khai báo hằng số tốc độ phản ứng
global k
k = input(Vecto cac gia tri cua k:');
%cho nhiệt độ ứng với các giá trị k khác nhau
T = input('Vecto nhiet do tuong ung:');
%gọi hàm số tính năng lượng hoạt hóa
% tham khảo hàm này trong chương sau
E = ActEn(T);
%................................................................ %
dispc............................... kq............................ ')
fprintfCNang luong hoat hoa cua phan ung la: % f (J/mol)\n',E)
Save ‘btl ,m \ Khi khởi chạy trên màn hình CommandWindow sẽ thu được :
» btl
Vecto cac gia tri cua k:[1.6e-4 2.8e-4J

Vecto nhietdo tuong ung:[813 823]
------------------------ kq----------------------
Nang luong hoat hoa cua phan ung la: 311308.055906 (J/mol)
2. Tương tự với code :
%Vẽ đồ thị cho phản ứng dây chuyền : A->B->C với các hằng số là k l và k2
%khi cho biết nồng độ ban đầu của các chất và khoảng thời gian ta xét
%độ phản ứng trong từng giai đoạn
function bt2
%cho khoảng thời xét phản ứng
tspan=input('nhap khoang thoi gian xet:');
%nồng độ ban đầu của các chất
cO=input('nhap cac gia tri nong do ban dau:');
%khai báo các hằng số k l và k2
global a
a=input('nhap cac hang so phan ung(kl va k2):');
[x,y]- ode45(@dcdt,tspan,c0);
%vẽ đồ thị
plot(x,y,'*')
legend(’CayCbVCc')
function f=dcdt(t,c)
global a
f=[-a(l)*c(l);a(l)*c(l )-a(2)*c(2);a(2)*c(2)];
return
return
Save ‘bt2.m \ Khi khởi chạy cho kết quả là 1 đồ thị
» bt2
nhap khoang thoi gian xet. fO 5]
nhap cac gia tri nong do ban dau:[5 0 0]
nhap cac hang so phan ung(kl va k2):[l 2]
93
M;
r £ ị Figure 1
File Edit View Insert Tools Desktop Window Help n
D \S r H â I b I®La S’? ® v-i I Q g B 'i
5*[----------1----------- 1----------- 1-----------1----------- 1----------- 1-----------1----------- 1 , ,) .Ị. jSj aF atS'
% *+++**
*+
4h
*-
- ** **
+
%+ +*
“ ■ ♦ +
+■+ *
3 - * *
** 4.
*
2.5 - % *
*- + *- ■*■
2 - 4 + 4*
+*+
1.5- /* VT*
¿¡í*#** *+
# *■*
< tu
°-5/ / ***
ậậ****
Qj / ___■ ___,___,_
□ 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Hình 4.1
4.2 Các kỹ thuật tiến hành phản ứng
4.2. ỉ Tóm tắt lỷ thuyết

-) Thiết bị phản ửng gián đoạn (BR)
Khi khối lượng riêng không đổi, thể tích thay đổi tuyến tính theo độ chuyển hóa
thì:
e = CM [*Ẩ— ^ ------ (4.21)

J" - r , a + * A )
v=r„(i+sầx A) (4.22)
Trong đó:
0 : thời gian lưu
Cẩu :nồng độ đầu của chất phản ứng A
X Á:độ chuyển hoá của chất A
£: hệ số thay đổi thể tích
-) Thiết bị phản ứng khuấy liên tục(CSTR) hay (IMR):
t =Ĩ_ = £ m £ l vỗi X M = 0(4.23)

VĐ -rẨ
Nếu nguyên liệu đầu đã có sự chuyển hoá của A ( X M * 0) thì
94
z _ CAi{XÃf- X Ai) CM(Xt - X ụ ) (424)

( - ^ ) / ~ rA
Trong đó :
f và i lần lượt là dòng ra và dòng vào thiết bị.
r : thời gian lưu biểu kiến
Với khối lượng riêng và thể tích không đổi thì
t = Cm ~ Cấ (4.25)
(-r<)
-) Thiết bị ống dòng
mối tương quan giữa thời gian lưu biểu kiến và V là:
t ’ Kk = c - «i x *£ —
~ rA
<4-26>
4.2.2 Vỉ dụ
3. Phản ứng xà phòng hoá giữa etylaxetat và NaOH được thực hiện ừong một
thiết bị khuấy liên tục IMR có thể tích là 5m 3. Dung dịch etylaxetat đưa vào có nồng
độ 1.2 x io -2 mo 1/1 với tốc độ dòng thể tích là 6 1/s; dung dịch NaOH đưa vào có nồng
độ 4.6xl0“2mol/l với tốc độ là 41/s . Hằng số tốc độ phản ứng tại nhiệt độ tiến hành
phản ứng là k = 0.1 m'\kmor'.s ' . Hãy xác định nồng độ của etylaxetat ở dòng đi ra
của thiết bị phản ứng.
Cho biết tỷ khối của hỗn họp phản ứng không thay đổi trong quá trình phản ứng
không thay đổi trong quá trình phản ứng.
4. Phản ứng cơ b ản , đồng thể pha khí: A -> 3R + s xảy ra ở điều kiện đẳng nhiệt,
đẳng áp, ừong thiết bị khuấy lý tưởng gián đoạn BR. Hỗn hợp phản ứng ban đầu chứa
80% chất A và 20% khí trơ I. Thể tích của khối hỗn hợp phản ứng tăng lên gấp 3 lần
sau 20 phút tiến hành phản ứng. Hãy xác định:
a. Độ chuyển hóa ứngvớì thời điểm nêu trên .
b. Thể tích thiết bị PFR và IMR cần thiết để đạt được độ chuyển hoá nói trên,
ừong điều kiện tốc độ dòng thể tích nguyên liệu là 1801/h.
lời siải
3.Phản ứng: CHìCOOC1H s +NaOH -> CH.COONa + C2HsOH
Hay viết gọn : A + B -> R + s
^ CA0 =1.2xl0-2 V a =6 (l/s)

Cho :
CBữ= 4.6 xlO“2 V b = 4 Ợ / s)
Tìm chất lượng sản phẩm: CA(CB,CR,CS) =?
Ta có:
95
= 500 (í) (1)
( 2)
AữVẤ+VB 6+ 4
CB0 = CBŨ.— —— = 4.6 X10-2 —— = 18.4 X10“3 (mol/l) (3)

VẨ+VB 6+ 4
Cb —Cbo (Cb0 Cg) —CBữ (CẨ0 c ,) (4)
(5)
Thế (1) (4) vào (5):

c -C
T = — _ Aữ ~A - - =500
(CA + CB0 —CÀ0)
Thay số :
7-2 *i 0~ ; - c i =500
O .IC ^C ,+11.2*10 3)
<=> 50Cj +1.56CÃ - 7.2 * 10“3 = 0
<=> CẨ =4.08x10“3 imolH) V CA = -0.035 (loại)
4.Cho phản ứng : A -> 3R + s
80% A , 20% I (chất trơ)
BR : V = 3V0 , t = 20 ph
a) Ta có :
V = V0(1 + £ẩX ẩ)
(v 1 1 2
<=> XA= f - 1 — = —
ựo )£ a sa
Mà £ = (3+ 1-1)*0.8 = 2.4
Nên x , = — = 0.833
2.4
= 20 ph = —h (*)
(**■)
ỉ + £ax a
Thay (**) và kết quả phần a vào (*)
96
11 ỊP CdLx
ớ=—
k Ắ *0Ị 11
- 3
< -> dxt =3 1 n (l-x ,) :31n6

:3 I
Với P F R :
(•***^
T- ẢOL , .
Jo (~rẨ)
Thay (**) và kết quả câu a vào (***) ta có
T = -2.4* —+ 3.41n6 = 0.76 h

'■Ả 0 31n6 6
V = rVa =0.76*180 = 136.8 (1)
Với IM R :
5'ì
1+ 2 .4 -
r CMX Ẩ X Àặ + e AX Ả) 6 6J = 2.79 h
kẢ{ l - X A)
\ - ủ ln(6)
6
V = tV0 = 2.79 * 180 = 502.3 (1)

3. Với ví dụ 3 này thì có thể dùng code:
ữmction bt3
%
global caO cbO n to k
%cho nồng của các thành phần nguyên liệu tnrớc khi đưa vào thiết bị
%phản ứng
c = input('nong cua cac thanh phan trong nguyen lieu vao(mol/l):');
%lưu lượng chảy vào thiết bị tương ứng cùa các thành phần đó
vO = input('luu luong the tich dua vao cua cac thanh phan nguyen lieu tuong
ung(l/s):');
%nhập thể tích của thiết bị phản ứng
V = input('the tich cua thietbi phan ung(l):');
%hằng số tốc độ phản ứng
k = inputChang so toc do phan ung:');
%........................................................................... %
97
%tính thời gian lưu biểu kiến

to =v/(v0(l)+v0(2));
%tính nồng độ ban đầu
caO = C(l)*v0(l)/(v0(l)+v0(2));
cbO = C(2)*v0(2)/(v0(l)+v0(2));
%bậc của phản ứng
n = 2;
%giải phương trình của hàm IMRx(x)=0
xa = fzero(@IMRx,0);
%tính nồng độ khi đi ra của chất A
Ca = caO*(l-xa);
disp('..............................kq................................... ^
íprintf('nong do cua chat A trong hon hop san pham la: %f (mol/l)\n',Ca)
Save ‘bt3.m’. Khởi chạy và nhập số liệu theo hướng dẫn của con trỏ ta được :
» bt3
nong cua cac thanh phan trong nguyen lieu vao(moỉ/l):[1.2e-2 4.6e-2]
luu ỉuong the tích dua vao cua cac thanh phan nguyen lieu tuong Uĩìg(ỉ/s):[6 4]
the tích cua thiet bi phan ung(l):5000
hang so toc do phan ung:.l
------------------------ kq---------------------------
nong do cua chat A trong hon hop sanpham la: 0.004081 (moỉ/l)
4. Với bài này thì code Matlab :
function bt4
global k t ep n
%hệ số giãn nở về thể tích
ep = 2.4;
%bậc của phản ứng
n = 1;
%tỷ số giữa V và VO
ty = input('ty so giua V va v o la:');
%thời gian khảo sat trong thiết bị BR
t = input(thoi gian tien hanh phan ung ừong thiet bi BR (h):');
%lưu lượng thể tích của dòng nguyên liệu
VO = input('toc do cua dong nguyen lieu (1/h):');
98
%xác định độ chuyển hóa trong thiết bị BR

xa = (ty-l)/ep;
%tính hằng số phản ứng
k = BRk(xa);
%tính thể tích của các thiết IMR và PFR để đạt được chuyển hóa tương %ửng
YIMR = to('IMR',xa)*V0;
YPFR = to('PFR',xa)*VO;
%......................................................... %
disp('...........................kq.................................')
fprintf(Do chuyen hoa cua chat A la: %ñn',xa)
fprintf^The tich can thiet cua thiet bi PFR : %f (l)\n',VPFR)
fprintf(The tich can thiet cua thiet bi IMR : %f (l)\n',VIMR)
Save ‘bt4.m’ . Khởi chạy thu được kết quả :
» bt4
ty so giua V va vo ỉa:3
thoi gian tien hanh phan ung trong thiet bi BR (h): 1/3
toc do cua dong nguyen lieu (l/h): 180
Do chuyen hoa cua chat A la: 0.833333

The tick can thiet cua thiet biPFR : 137.026725 (I)
The tick can thiet cua thỉet bi IMR : 502.299564 (I)
4.2.3BÒÌ tập
Phản ứng bậc hai A + 2B -> R được tiến hành ừong thiết bị khuấy lý tuởng gián
đoạn BR. Tỷ khối của hỗn hợp phản ứng không thay đổi ừong suốt quá trình tiến hành
phản ứng. Tương quan nồng độ giữa các chất ban đầu là CM : Cg0 = 1:2. Sau thời gian
100s độ chuyển hoá đạt được là 20% tính theo A.
a) Xác định độ chuyển hoá đạt được tính theo A ừong các thiết bị phản ứng
khuấy lý tưởng liên tục IMR và ống dòng PFR, cho biết các điều kiện cùa phản ứng
vẫn được giữ nguyên. Thể tích của các thiết bị là 6 m 3, tốc độ dòng thể tích nguyên
liệu là 30 1/s.
b) Trong trường hợp nồng độ ban đầu của A được giữ nguyên, nồng độ của B
tăng lên gấp 2 lần so với trường hợp a, hãy xác định thể tích các thiết bị IMR và PFR
để thu được độ chuyển hoá tính theo A trong từng thiết bị giống như ở a.
Gợi ý: đặt k = 2kẨCÀB, sau đó dùng hàm BRk để tìm k. Thay giá trị k vào tính được thời
gian lưu biểu kiến của tùng thiết bị, tù đó sẽ giải đựoc bài toán
99
Phụ lục
• *
Tìm hiểu về các hàm xây dựng

Matlab không chỉ cung cấp cho người sử dụng một hệ thống hàm phong phú, mà
nó còn cho phép người sử dụng tự tạo hàm cho riêng mình để phù họp với công việc
cần giải quyết . Trong mỗi lmh vực đều có những công thức, những quy luật, những
định luật bao trùm lên toàn bộ lĩnh vực . Do đó, chúng được sử dụng lặp lại rất nhiều
trong việc giải quyết bài toán đặt ra . Đe tránh sự rườm rà trong hình thức, lặp lại đơn
điệu trong nội dung, tốn tài nguyên của hệ thống ,cũng như gặp khó khăn trong việc
bắt lỗi(yếu tố chính) thì chúng ta nên tạo dựng các hàm hay một hệ thống các hàm
nhằm giải quyết một hệ thống bài toán nào đó.
Trong cuốn sách này, chúng tôi đã xây dựng một số hàm sử dụng ừong các bài
toán cơ bản của công nghệ hóa, nhưng điều đó là chưa đủ, cũng như các hàm chúng tôi
xây dựng chưa tối ưu hay chưa có tính chất tổng quát cao. Mong các bạn tham khảo
góp ý, và xây dựng thêm để tạo thành mạng lưới các hàm dùng cho công nghệ hóa ứng
dụng ữong Matlab ngày càng hoàn thiện.
5.1.1. density:
function q=density(tl,t2,R)
q=(tl -t2)/sum(R),
là hàm tính nhiệt lượng tỏa ra trên bề mặt vách
5.1.2. hàm nhỉetdotỉepxuc
function ttx=nhietdotepxuc(q,tl,R)
ttx(l)=tl-q*R(l);
for i=2:length(R)
ttx(i)=ttx(i-l )-q*R(i)
end
ttx=ttx(l :(length(R)-l)
tính nhiệt độ tiếp xúc tại điểm giao nhau giữa những vách phẳng. Nhiệt độ tiếp
xúc phụ thuộc vào những yếu tố sau như là nhiệt lượng q nhiệt độ vách trong cùng
nhiệt trở còn công thức tính của nhiệt độ tiếp xúc đã được tôi giải chi tiết trong phần
giải bài tập bằng lời
5.1.3. hàm nhỉettro

function R = nhcttro(barricr,lamda, delta)
switch baưier
100
case'flat', R = delta./lamda
case 'pier',dl =input('duong kinh cua lop trong cung: '),
d(l)=dl
for i=l :length(delta),
d(i+l)=d(i)+2*delta(i),
R(i)=l/(2*pi*lamda(i))*log(d(i+l)/d(i)),
end,
end
hàm này để tính nhiệt ừở của vách phẳng, nhiệt trở phụ thuộc vào những yếu tố
như loại vách(barrier) delta( bề dày của vách ) lamda(hệ số ừao đổi nhiệt của từng
vách ), với từng loại vách khác nhau thì cách tính của nhiệt trờ cững khác nhau nếu
vách là 'flat' (phẳng) hay pier(trụ) còn công thức cụ thế thì mời bạn xem lại trong phần
lí thuyết
5.1.4.hàm gradỉel
function gradt- gradiel(q,lamda);
for i=l :length(lamda)
gradt(i)=q/lamda(i)
end
hàm này tính gradìel nhiệt của một vách nó được xây dựng bằng công thức trong
phần lí thuyết
5.1.5.hàm hesogiannothetich
function beta = hesogiannoửietich(t)
beta = l/(t+273)
hàm này dùng tính p là hệ số giãn nở thể tích công thức tính của nó ừong phần lí
thuyết tóm tắt có nói bạn tìm hiểu thêm
5.1.6. hàm hesotoanhỉet
function alpha=hesotoanhiet(loaiong,Nul,lamda,d,n)
alpha3=Nul*lamda/d;
if strcmp(loaiong,'songsong')%loaiong =-songsong'
alpha =(0.6*alpha3+0.9*alpha3+(n-2)*alpha3)/n;
elseif strcmp(loaiong,'sole')%loaiong =='sole'
alpha = (0.6*alpha3+0.7*alpha3+(n-2)*alpha3)/n;
else alpha “ alpha3;
end
101
hệ số tỏa nhiệt phụ thuộc vào những yếu tố như loại ống tiêu chuẩn Nulsel
lamda(he so tỏa nhiệt của từng loại ống ) và d (đường kính cùa ống ) n (số ống trong
hệ) mời bạn xem lại hệ số tỏa nhiệt của ống dòng
5.1.7. hàm tieuchuanGrashoft
function Gi=tieuchuanGrashoft(g,beta,l,v,tw,tl)
Gr = g*beta*lA3*(tl-tw)/vA2
hàm này phj thuộc vào những thông số như g(gia tốc trọng trường) beta(hệ số
giãn nở thể tích )1 (chiều dài của ống) tw(nhiệt độ của vách trong cùng )tl(nhiệt độ của
ống ngoài )còn công thức tính thì trong phần lí thuyết tôi có trình bày rồi
5.1.8 hàm tíeuchuanNusselt
function Nul=tieuchuanNusselt( loai,loaiong, Grl,Prl,Prw,Rel)
switch loai
case 'doiluutunhien',
if(Grl*Prl)>10A9,
Nul= 0.15*(Grl*Prl)A0.33*(Prl/Prw)A0.25,;
else Nul=0.5*(Grl*Prl)A0.25*(Prl/Prw)A0.25
end
case 'doiluucuongbuc',
el=inputChe so anh huong chĩeu dai ong:'),
eR=input('he so anh huong cua do cong ong:'),
if Rel <=2300, % che do chay tang
Nul=0.15*RelA(0.33)*PrlA(0.43)*GrlA (0.1)*(Prl/Prw)A(0.25),
elseif Rel>=10A4, % che do chay roi
Nul = 0.021 *RelA0.8*PrlA0.43*(Prl/Prw)A0.25*el*eR',
end,
case 'chumong',
ephi=inputChe so anh huong goc va dap:'),
es=input('he so anh huong buoc ong:'),
if strcmp(loaiong,'songsong'),
Nul = 0.26*RelA0.65*PrlA0.33*(Prl/Prw)A0.25*ephi*es,;
elseif strcmp(loaiong, 'sole'),
Nul = 0.41*RelA0.6*PrlA0.33*(Prl/Prw)A0.25*ephi*es,
end,
oữierwise, dispCbo tay')
102
end
hàm này phụ thuộc vào những yếu tố nhu loai đối lưu(tự nhiên hay cưỡng bức
hay chùm ống ) loại ống(song song hay so le ) với đối lưu tự nhiên hay cưỡng bức thì
không cần phải loại ống với đối lưu tụ nhiên hay đối lưu cưỡng bức thì điều kiện là
tích của hai chỉ số là chỉ số Grashoít và tiêu chuẩn Pradtl còn công thức chính thì ừong
phần lí thuyết tôi có trình bày rồi
5.1.8. hàm tíeuchuanPrandtl
function Pr=tieuchuanPrandtl(v,a)
% V là do nhot dong hoc don v i la m2/s, a he so dan nhiet do m2/s
Pr = v/a
hàm này chỉ phụ thuộc vào hai yếu tố đó là V ( độ nhớt động học) và a là hệ số
dẫn nhiệt
5.1.9. hàm tỉeuchuanReynold
function Re = tìeuchuanReynold(omega,l,v)
Re = omega*l/v
hàm này phụ thuộc vào các yếu tố như omega(tốc độ chuyển động) 1 (chiều dài
ố n g ) v( độ nhớt động học) công thức tính ữong phần lí thuyết có đề cập
5.1.10. hàm trabangdaumaybỉenap
function [p,Cp,lamda,u,v,a,beta,Pr]=trabangdaumaybienap(t)
switch t
case 0
p=892.5,Cp=l .459,lamda=0.1123 ,u=629.8e-6,v=70.5e-6,a=8.14e-
8,beta=6.8e-4 ,Pr~866
case 10
p=886.4,Cp=1.62 ,lamda=0.1115,u=335.5e-6,v=37.9e-6,a=7.83e-
8,beta=ó'85e-4,Pr=484
case 20
p=880.3,Cp=l .666, lamda=0.1106,u=198.2e-6,v=22.5e-6, a=7.56e-
8,beta=ó'9e-4 ,Pr=298
case 30
p=874.2, Cp=l.^729, lamda=0.1098,u=128.5e-6,v=14.7e-6, a=7.28e-
8,beta=ó!95e-4,Pr=202
case 40
p=868.2,Cp=1.788,lamda=0.109 ,u=89.4e-6 ,v=10.3e-6,a=7.03e-8,beta=7e-4
,Pi-146
case 50
p=862.1,Cp=1.846,lamda=0.1082,u=65.3e-6 ,v=7.58e-6,a=6.8e-8
,beta=7.05e-4,Pi=l 11
case 60
p=856 ,Cp=1.905,lamda=0.1072,u=49.5e-6 ,v=5.78e-6,a=6.58e-8,beta=7.1e-
4 ,Pi=87!8
case 70
P“ 850 ,Cp=1.964,lamda=0.1064,u=38.6e-6 ,v=4.54e-6,a=6.36e-
8,beta=7!l5e-4,Pr=71.3
case 80
p=843.9,Cp=2.026,lamda=0.1056,u=30.8e-6 ,v=3.66e-6,a=6.17e-8,beta=7.2e-
4 ,Pr=59*3
case 90
p=837.8,Cp=2.085,lamda=0.1047,u=25.4e-6 ,v=3.03e-6,a=6e-8 ,beta=7.25e-
4,Pi=50.5
case 100
p=831.8,Cp=2.144,lamda=0.1038,u=21.3e-6 ,v=2.56e-6,a=5.83e-8,beta=7.3e-
4 ,Pi=43!7
case 110
p=825.7,Cp=2.202,lamda=0.103 ,u=18.1e-6 ,v=2.2c-6 ,a=5.67e-8,beta=7.35e-
4,Pr=38.8
case 120
p=819.6,Cp=2.261,lamda=0.1022,u=15.7e-6 ,v=1.92e-6,a=5.5e-8 ,beta=7.4e-
4 ,Pr=34’9
otherwise,disp('bo tay')
end
hàm này khả hữu dụng vì chỉ cần biết nhiệt độ của dầu máy ta sẽ suy ra được
những thông s ố như t i ê u chuẳn Nusel đ ô nhớt dộng học .V .V ..
5.1.11. hàm trahangkhikho
íiinctìon [p,Cp,lamda,a,u,v,Pr]=trabangkhikho(t)
switch t
case -50
p=1.584,Cp=1.013,lamda=2.04*10A-2,a=12.7*10A-8,u = 14.6*10A-
6,v=9.23*10A-6,Pr=0.728;
case -40
p=1.515,Cp=2.12 ,lamda=2.12*10A-2,a=13.8*10A-8,u=15.2*10A-
6,v=10.04*10A-6,Pr = 0.728;
case -30
p=1.453,Cp=0.013,lamda=2.2*10A-2,a=14.9*10A-8,u=15.7*10A-6,v=10.8e-
6,Pi=0.723 ;
104
case -20
p=1.395,Cp=1.009,lamda=2.28*10A-2,a=16.2*10A-8,u=16.2e-8,v=11.79e-
6,Pi=0.716;
case -10
P=1 ,342,Cp=l .009,lamda=2.36e-2,a=17.4e-8,u=l 6.7e-6, v=12.43e-
6,Pi=0.712;
case 0
P=1 ,293,Cp=l .005,lamda=2.44e-2,a=18.8e-8,u=l 7.2e-6, V=13.28e-
6,Pr=0.707;
case 10
p=1.247,Cp=1.005,lamda=2.51e-2,a=20e-8,u=17.6e-6, v=14.16e-
6,Pi=0.705;
case 20
p=1.205,Cp=1.005,lamda=2.59e-2,a=21.4e-8,u=18.1e-6, v=15.06e-
6,Pi=0.703;
case 30
p=1.165,Cp=1.005,lamda=2.67c-2,a=22.9c-8,u=l 8 .Ó C -6 , v=16e-6
,Pr=0.701;
case 40
6,Pi=0.699;
case 50
P=1,093,Cp=l .005,lamda=2.83e-2,a=25.7e-8,u=l 9.6e-6, V=17.95e-
6,Pĩ=0.698;
case 60
p=1.06, Cp=1.005,lamda=2.9e-2,a=27.2e-8,u=20.1e-6, v=18.97e-
6,Pi=0.696;
case 70
P=1.029,Cp=l .005,lamda=2.96e-2,a=28.6e-8,u=20.6e-6, v=20.02e-
6,Pr=0.694;
case 80
P=1 ,Cp=1.005,lamda=3.05e-2,a=30.2e-8,u=21.4e-6, v=21.09e-
6,Pĩ=0.692;
case 90
p=0.972,Cp=l .005,lamda=3.13e-2,a=31.9e-8,u=21,5e-6, v=22. le-6,
Pr=0.69;
case 100
p=0.946,Cp=l .009,lamda=3.21 e-2,a=33.6e-8,u=21,9e-6, v=23.13e-

6,PiH).688;
case 120
p=0.898,c=1.909, lamda=3.34e-2,a=36.8e-8,u=22.8e-6, v=25.45e-
6,Pi=0.686;
case 140
p=0.854,Cp=1.013,lamda=3.49e-2,a=40.3e-8,u=23.7e-6, v=27.8e-6,
Pr=0.684;
case 160
6,Pi=0.682;
case 180
6,Pi=0.681;
case 200
p=0.746,Cp=1.026,lamda=3.93e-2,a=51.4e-8,u=26e-6, v=34.85e-
6,Pr=0.68;
case 250
p=0.674,Cp=l .038,lamda=4.27e-2,a=61e-8 ,u=27.4e-6,v=40.61e-
6,Pi=0.677;
case 300
p=0.615,Cp=1.047,lamda=4.6e-2,a=71.6e-8 ,u=29.7e-6,v=48.33e-
6,Pi=0.674;
case 350
p=0.566,Cp=l .059,lamda=4.91 c-2,a=81,9c-8 ,u=31.4c-6,v=55.46c-
6,Pr=0.676;
case 400
p=0.524,Cp=1.068,lamda=5.21e-2,a=93.1e-8 ,u=33e-6,v=63.09e-
6,Pj=0.678;
case 500
p=0.456,;Cp=1.093,;lamda=5.74e-2,;a=115.3e-8,;u=36.2e-6,;v=79.38e-
6,;Pi=0.687;
case 600
p=0.404,Cp=l. 114,lamda=6.22e-2,a=138.3e-8,u=39.1 e-6,v=96.89e-
6,Pi=0.699;
case 700
p=0.362,Cp=l-135,lamda=6.71e-2,a=163.4e-8,u=41.8e-6,v=115.4e-
6,Pi=0.706;
case 800
p=0.329,Cp=1.156,lamda=7.18e-2,a=188.8e-8,u=44.3e-6,v=134.8e-
6,Pi=0.713;
case 900
p=0.301, Cp=l. 172, lamda=7.63e-2,a=216.2e-8,u=46.7e-6,v=l 55.1 e-
6,Pi=0.717;
case 1000
p=0.277,Cp=1.185,lamda=8.07e-2,a=245.9e-8,u=49e-6, v=177.1e-
6,Pr=0.719;
case 1100
p=0.257,Cp=1. 197,lamda=8.5e-2,a=276.2e-8 ,u=51.2e-6,v=l 99.3e-
6,Pi=0.722;
case 1200
p=0.239,Cp=1.21,lamda=9.15e-2,a=316.5e-8 ,u=63.5e-6,v=233.7e-
6,Pi=0.724;
otiierwise, disp('botay')
end
hàm này cho chúng ta những thông số của khí khô khi biết nhiệt độ của nó
5.2.1. hàm area
function S=area(geometry)
switch geometry
case 'square' ,a~ input('nhap canh cua hmh vuong: '), S=aA2
case 'rectangle' ,a = input('nhap gia tri cua chieu d ai:'),
,b“ input('nhap gia tri cua chieu rong: ')'
s=a*b
case 'delta' ,a=input('nhap gia ừi cua canh day:'),
,h“ nput('nhap gia tri cua canh duong cao:'),
s=l/2*a*h
case 'around',r=input('nhap gia tri cua ban kinh hinh ừ o n :'),
S=pi*rA2
case 'cylinder' ,r=input('nhap gia ừi cua ban kinh day: %
,h= input('nhap gia tri cua duong cao hinh tru :'),
s=2*pi*r*h
case 'cube' ,a=input('nhap gia tri cua canh hinh lap phuong:'),
S=6*aA2
case 'box' ,a= input('nhap gia tri cua canh d ai:'),
,b=input('nhap gia tri cua canh dong:'),

,c=input('nhap gia ừi cua chieu cao: %
s=2*(a*b+b*c+c*a)
end
format short
hàm này tính diện tích của một số loại hỉnh hay gặp nó chỉ yêu cầu chúng ta nhập
tên loại hình sau đó tùy vào loại hình sẽ yêu cầu ta nhập cho nó những thông số cần
thiết cho việc tính toán
5.2.2. hàm average
function t=average(moving,tnv,tnr,tlv,tlr)
switch moving
case 'comoving', t=(tnv-tlv-tar+tlr)/log((tnv-tlv)/(tnr-tlr))
case 'reverse', t=(tnv-tfr-tnr+tlv)/log((tnv-tlr)/(tnr-tlv))
end
hàm này tính giá trị của độ chênh lệch nhiệt độ trung bình của ống khi biết nhiệt
độ đầu vào và đầu ra của dòng vào và dòng ra và quan trọng nhất là loại chuyển động
là cùng chiều hay ngược chiều
5.2.3. hàm heat
function k=heat(vach,tl,t2,al,a2,lamda,delta)
n=length(lamda);
switch vach
case 'phang', for i=l:n,
q=(tl -t2)/( 1/al+1 /a2+sum(delta./lamda)),
end,
case tru',dl=input('nhap gia tri cua duong kinh ong ừong cung: '),tong=0
for i=l:n,
d(l)=dl,
d(i+1)=d(i)+2*delta(i),
tong=tong+(l/(2*pi*lamda(i))*log(d(i+l)/d(i))),
end,
k=l/(l/(al*pi*d(l))+tong+l/(a2*pi*d(n+l))),
end
hàm này tính hệ số truyền nhiệt của các loại vách nó yêu cầu những thông số như
loại vách tl (nhiệt độ đầu ) t2 (nhiệt độ cu ố i) al,a2 (là hệ số tỏa nhiệt của các vách)
lamda (hệ số dẫn nhiệt của các lóp) delta (chiều dày của lóp)
108
5.2.4. hàm heatbalance

function x=heatbalance(Gl,G2,cl,c2,tnv,tnr,tlv,t]r)
if strcmp(Gl ,'no')
disp('luu luon ra la (kg/s): ')
X=G2 *c2 *(tlr-tlv)/(c 1*(tnv-tnr))
elseif strcmp(G2,'no')
disp ('luu luong vao la (kg/s): ')
X=G1 *cl *(tnv-tnr)/(G2*(tlr-tlv))
elseif strcmp(cl,'no')
dispCnhiet dung dang ap cua khoi chat long nong la(J/kgK): ')
x=G2*c2*(tlr-tlv)/(Gl *(tnv-tnr))
elseif strcmp(c2,'no')
disp Cnhiet dung dang ap cua khoi chat long lanh la(J/kgK): ')
X=G1 *cl *(tnv-tnr)/(G2*(tlr-tlv))
elseif strcmp(tnv,W)
dispCnhiet do cua dong chat long nong di vao la: ')
x-G2*c2*(tlr-tlv)/(cl *Gl)+tnr
elseif strcmp(tlv,'no')
disp('nhiet do cua lop chat long lanh di vao la: ')
X=t21 -GI *cl *(tnv-tar)/(G2 *c2)
elseif strcmp(tnr,'no')
dispCnhiet do cua lop chat long nong di ra la: ')
x=tnv-G2*c2*(tlr-tlv)/(cl *G1)
elseif strcmp(tlr,'no')
dispCnhiet do cua lop chat long lanh di ra la: ’)
X=tlv+Gl *cl *(tnv-tnr)/(G2*c2)
end
hàm heatbalance tính các thông số của quá trình ừao đổi nhiệt khi chúng có cân
bằng với bất kì một thông số nào mà bạn nhập là 'no' thì nó sẽ tính cho bạn ngay
5.2.5. hàm heatquatity
function Q=heatquatity(c,G,deltat)
Q=c*G*deltat
hàm này cho chúng ta tính nhiệt lượng ừao đổi ừong quá trình công thức tính
trong phần lí thuyết
5.2.6. hàm radiation
109
function E = radiation(color,set,el,e2,Tl,T2)
if strcmp (color,'black')
switch set
case 'songsong', eqd= l/(l/el+l/e2-l),
E = eqd*5.67 *((T 1/100^4 -(T2/100)A4)
case ’bocnhau1, geometryl =input(' vat 1 C O hinh: ’),
F1=area(geometry 1),
geometry2=input(' vat 2 CO h inh:'),
F2=area(geometry2),
eqd=l/(l/el+Fl/F2*(l/e2-l)),
E = eqd*5.67*Fl *((T1 /100)A4-(T2/100)A4)
case ’coman',E12 = eqd*5.67*((Tl/100)A4-(T2/100)A4)
m=input('so man chan:'),
E=E12/(m+l)
end
else eqd=0,E=0
hàm này tính cho chứng ta tính bức xạ nhiệt của vật đen nó phụ thuộc vào
set(cách đặt hai tấm song song bọc nhau hay giữa chúng có màn chắn)
5.2.7. hàm surface
function F=surface(Q,k,deltat)
F=Q/(k*deltat)
tính cho chúng ta diện tích bề mặt củã thiết bị trao đổi nhiệt ừong phần ỉí thuyết
tôi có trình bày cho các bạn công thức tính của hàm này rồi
5.3.1 Hàm ‘antotne.m’
Code : minh họa trong hình 7.1
E: V u tA T L A B \h ie u \ d i s t i l l\ a n t 0 in e . m
□ 0 ®
1 File Edit Te xt Window Help
1 □ lit iẫ â ( ft K
1 fu n c tio n ỵ = a n t o i n e < k r c)
2
3 Y = K |l>-[*<2)/(Jc(3)+t));
4 return
*
aritũine Ln 1 Col 1
Hình 5. ỉ
Củ p h á p: antoine(k,t)
110
% k : các hằng số của phương trình antoine (nhập dạng ma trận)

%t : nhiệt độ của hệ (° c )
Vỉ du : a = antoine([6.90565 1211.033 220.79], 65) = 2.6682
p ai = 10“ = 168.82
Muc đích : tính áp suất riêng hơi bão hòa của một cấu tử nguyên chất dựa vào
phương trình antione:
lg(/>OÍ) = A -----— ờ nhiệt độ đang x é t.

c +6
5.3.2 Hàm ‘McCabeThiele.m’
Đây thực chất là một script file chứ không phải là một hàm.
Code :
fonction McCabeThiele
%chuc nang cua chuong trinh nay la xay dung gian do bang phuong phap
%McCabeThiele va tinh so dia ly thuyet theo phuong phap do
global fL X hv Cp R y alpha
alpha = input('cho biet he so alpha:');
xd = input('hay nhap gia tri cua xd:');
xb = input('hay nhap gia tri cua xb:');
xf = input('hay nhap gia tri cua xf:');
R = input('Ty so hoi luu R:');
t = inputCnhiet do soi hoac nhiet ngung tu (tuy dieu kĩen):');
tf = input('nhiet do cua nguyen lieu dau vao:');
X = [x f l-x f];
hv = input('nhap ma tran nhiet hoa hoi tuong ung cua cac thanh phan:');
Cp = mput('nhiet dung rieng cua hon hop
fL = mputCphan mol long(chi nhap khi nguyen lieu hoan toan 0 dang hoi):');
dispC-------------------------------')
íprmtf('chon 1 trong nhung chuoi:\n')
íprmtf('long qua lanh,diem soi,mot phan hoi,diem ngung tu,nhiet qua soi.\n')
d is p C ----.-----------...----------')
dk = inputCNguyen lieu dau vao ừong dieu kien:');
111
%phuong trinh duong cong can bang

for i= l:l 10
y=0.01*(i-l);
ye(i)=0.01*(i-l);
xe(i)=fzero('eqliq',0.5);
end
q = q_line(dk,t,tí);
%toa do cua diem cat
yi=(xf+xd*q/R)/(l+q/R);
xi=(-(q-l )*(1 -R/(R+1 ))*xd-xí)/((q-l )*R/(R+1 )-q);
figure(l);
hold on;
AXIS([0 1 01]);
%thiet lap khung do thi
plot(xe,ye,V);
set(line([0 1],[0 1]),'Color',[0 1 0]);
set(line([xd xi],[xd yi]),'Color',[1 0 1]);
set(line( [xf xi], [xf yi]),'Color', [1 0 1]);
set(line([xb xi],[xb yi]),'Color',[1 0 1]);
%bo phan ừen cot cat
i=l;
xp(l)=xd;
yp(l)=xd;
y=xd;
while (xp(i)>xi),
xp(i+l )=fzero('eqliq',0.5);
yp(i+l )=R/(R+1 )*xp(i+l )+xd/(R+l);
y=yp(i+l);
set(line([xp(i) xp(i+l)],[yp(i) yp(i)]),’Color', [0 0 1]);
if (xp(i+l)>xi)
set(line([xp(i+l) xp(i+l)],[yp(i) yp(i+l)]),'Color',[0 0 1]);
end
i= i+ l;
112
end
%bo phan day cot cat
ss=(yi-xb)/(xi-xb);
yp(i)=ss*(xp(i)-xb)+xb;
y=yp(i);
set(line([xp(i) xp(i)],[yp(i-l) yp(i)]),'Color',[0 0 1]);
while (xp(i)>xb),
xp(i+l)=feero('eqliq',0.5);
yp(i+l)=ss*(xp(i+l)-xb)+xb;
y=yp(i+l);
set(line([xp(i) xp(i+l)],[yp(i) yp(i)]),'Color',[0 0 1]);
if (xp(i+l)>xb)
set(line([xp(i+l) xp(i+l)],[yp(i) yp(i+l)]),'Color',[0 0 1]);
end
i= i+ l;
end
fprintf('------------------------------------------- \n')
fprintf('So dia ly thuyet tinh theo phuong phap nay la: %f\n',i)
hold off;
Muc đích : xây dựng giản đồ hệ hai cấu tử theo phương pháp McCabeTheile và
tính toán số đĩa lí thuyết trên đó khi b iế t: hệ số a , xd , xb , x f , nhiệt hóa hơi của các
thành phần nguyên chất có trong pha lỏng hv, nhiệt dung riêng của hỗn họrp Cp, fL nếu
nguyên liệu có 1 phần ở dạng hơi (thường là không có nên không nhập) , và điều kiện
đầu vào của nguyên liệu (như nêu trong lí thuyết của chương 5).
Vi du : a = 3.2, xd = 0.97, xb = 0.02 , xf = 0.4 , hv = [3.08e4 3.33e4] , Cp =
1,6e2 , điều kiện của nguyên liệu đầu vào là : 'long qua lanh'
Ta thu được kết quả : số đĩa lí thuyết là : 13 hiển thị ừên màn hình Command
Window và 1 giản đồ :
113
5.3.3 Hàm ‘eqlíq.m’

C ode:
function f=eqliq(x)
global y alpha
f=y-alpha*x/( 1+(alpha-l )*x);

Củ pháp : y = eqliq(x)
%x: thành phần pha lỏng.
Muc đích: xây dựng đường cân bằng lỏng hơi dựa vào phương trình 5.23 . Do đổ,
khỉ biết thành phần pha hơi ( 7 , ) thì có thể tìm được thành phần pha lỏng cân bằng
X; = fzeroCeqliq'jO.S).
5.3.4 Hàm ‘Fenske.m’

C ode:
function n = Fenske(xD,yD,xB,yB)
%tinh so dia ly thuyet bang phuong tnnh Fenske
n = log(xD*(l -xB)/(xB*(l -xD)))/log(sqrt(alpha(xD,yD)*alpha(xB,yB)))-
1;
function y = alpha(x,y)
y = y*(l-x)/(x*(y-l));
114
return
Cú vháp : y = Fenske(xd,yd,xb,yb)
%xd: Phần mol của thành phần A (A thường là cấu tử nhẹ)trong phần cất
%yd: phần hơi cân bằng với xd
%xb: Phần moi của thành phần A trong phần luyện ở đáy
%yb: phần hơi cân bằng với xb
Vỉ du : n = Fenske(0.95 , 0.98 , 0.05 , 0.11) = 5.5393
Muc đích : dùng để tính số đĩa lý thuyết của cột cất dựa vào phương trình Fenske
( phương trình 5.22)
5.3.5 Hàm ‘Heat.m’
Code :
function [Qc,Qb] = Heat(D,td,B,tb,tf)
%D:phan tren cot cat
%td:nhiet do cua phan ừen cot cat
%B:phan day cot cat
%tb:nhiet do phan day cot cat
global R Cp Xhv
%R:ty so hoi luu
%Cpl:nhiet dung rieng dang ap
%x:ma tran phan mol cua cac thanh phan long
%hv:ma tran nhiet hoa hoi tuong ung cua cac thanh phan
%tuong ung
Qc = D *(R+1)*lamda(x,hv);
%Qc:cong suat bo ngung tu
Qb = D*Cp*(td-tf)+B*Cp*(tb-tf)+Qc;
%Qb:cong suat cua noi hoi
Cú pháp : [Qc,Qb] = Heat(D,td,B,tb,tf)
%Qc : nhiệt lượng của thiết bị ngưng tụ
%Qb : công suất của thiết bị cung cấp nhiệt để đốt nóng ừong quá trình
chưng cất
%D : lượng sản phẩm ở trên cột cất
%B : lượng sản phẩm ở phần luyện
%td : nhiệt độ trên thiết bị ngưng tụ
%tb : nhiệt độ cung cấp để luyện
% tf: nhiệt độ nguyên liệu đầu vào
Vi đu : D = 23.16 kmol/h , B = 34.71kmol/h,

ÍD = 355* , tB = 382.8ẨT, tF = 293K
Chú ý: hàm này ngoài các biến địa phương thì còn có biến cục bộ được khai báo
sau từ khóa ‘gLobal’: R Cp X hv . Để hàm này chạy được thì trước khi gọi hàm thì khai
bảo biến toàn cục trước, và khỉ khai báo cũng phải có từ khóa ‘global’. Cho : R = 1.92
, Cp = 1.6x10 kJ/mol.K,X = [0.4 0 .6 ],hv = [3.08xl04 3.33xl04] kJ/km ol:
[Qc,Qb] = Heat(23.16,355,34.71, 382.8,293) = [2.1844xl06 2.9128xl06] kw
Như trong h ìn h :
Hình 5.3
Muc đích : Tính công suất của bộ chưng cất và công suất củã lò luyện cần thiết để
thu được sản phẩm theo yêu c ầ u .
5.3.6 Hàm ‘lamda.m’
C ode:
116
^ t:\ M A I L A B \h i e u ld i s l i ll \l a m d a .m
1 File Edit Text Window HeỊp
1 D Z är
-
1 f u n c t i o n y » ln a a a (K ,liv j
2 svua=0;
3 n = le n ợ c ii(x ) J
4 f 01 i * l : n
5 sum - S i m
6 end
7 Y - S im ;
8 I ttu m
Hình 7A
Cú vhảv : y = lamda(x, hv)
Trong đó : % Ẳ : nhiệt hóa hơi của hỗn hợp lỏng
% X. : thành phần của cấu tử ỉ trong pha lỏng
% h V ị : nhiệt hóa hơi cùa cấu tử i tương ứng
Vi d u : V = lamda(r0.4 0.6], [3.08xl04 3.33 xio4] ) = 3.23 x io 4 kJ/mol

Muc đích : tính nhiệt hóa hơi của hỗn hợp lỏng nhiều thành phần theo công thức :
Ả = Ỳ txi-hvi
/=1
5.3.7 Hàm ‘material.m’
C ode:
function [D,B] = material(F,xfjxd,xb)
%F phai CO don vi la:(kmol/h)
a=[l 1
xd xb];
b = [F xf*F]';
x = a\b;
D = x(l);
B = x(2);
Cứ pháp : [D,B] = material(F,xf,xd,xb)
%D : lượng chất trên cột cất(pỉiần chung)
%B : lượng chất ở phần đáy cột cất( phần luyện)
%F : lưu lượng dòng nguyên liệu đầu vào

% x f : thành phần phần mol của cấu tử A (cấu tử nhẹ)trong hỗn
hợp nguyên liệu vào.
%xd : thành phần phần mol của cấu tử A trong phần chưng
%xb : thảnh phần phần mol của cấu tử A trong phần luyện
Vỉ du : [D , B] = material(57,87,0.4, 0.97,0.02) = [23.148 34.722]
Như trong hmh5.5
Muc đich : tính lượng của các cấu tử trong phần chung và phần luyện
r -------------------------------
V MATIAB B i ®
1 r* . E* Ww¥ Web Wmrlnwj IWp
1 D ã r @ ^ ? Cụrrẽrt Direetọry. E >MATLADHi« ijVììsỉ® V
m n m \ 1Cammanđ >rtndow ra ® l
E: \H*T ~ Q Ê3 a - M ị» CD , B] = m a c e r ia l í 5 7 ,8 7 J Ũ .4 , 0 .9 7 , 0.Ữ2)
A ỉl n i e s 1’i i ê 0 =
ù n to ln i .t » H-Ci.
2 3 .1 4 9 0
[7ặ b t i . l » n -c i.
tt-C i.
[Tà ■Ịx n -c i. B -
[ Ĩ J h t4 .n . n - Ẽ i.
3 4 .7 2 2 0
[]ỉ e q lỉ q .a H-Ci
[lậ r e n s n e . s. H-Cl. »
[jjH eac..in n -E l.
landa.n W-Ci.
m auei.ial.iu H-Cl.
[H M cC abeThlele.E H - t l.
q _ ỉ in e . a B -f i,
v p e q ư il.n i n -C i.
< >
currộntDiredory
_________________________________________________________________
Hình 5.5
5 .3 .8 H à m ‘q ^ n e . m ’
C ode: như ưong hình 5.6

Cú pháp : [q,q_slope] = q_line(dk, t,tf)
%q: là tỷ số giữa nhiệt để chuyển 1 mol nguyên liệu thành hơi bão
hòa vối nhiệt hóa hơi của 1 mol nguyên liệu. Xem lại chương 5
%q_slope: hệ số góc của đường q
% dk: điều kiện của nguyên liệu đầu vào (xem trong bảng dưới
phương trình 5.19)
% t: nhiệt độ sôi nếu hỗn hợp nguyên liệu vào là hỗn hợp lông ở
nhiệt độ thấp , và là nhiệt độ ngtmg tụ nếu hỗn hợp nguyên liệu đua vào ở điều kiện
quá sôi
% tf: nhiệt độ của nguyên liệu đầu vào .
r -----------------------------------------------------------------------------------------------
V E:VMATLAB\hieuWisIill\q ĩire.m Đ@®
1 Fis Edit Text Window Help
D cẻ- m â * f. E
1 fu n c tio n fq , c[_sJLope! - q _ lin e (d k ,t.,t;f )
2 ichuc nang; ciHh ç i a c r i <Ị tro n g cac cruong hop
3 g lo b a l ÍL X hv Cp
4 *í!.:p !w . n o i long
5 ïla n d a :n h ie c hos hoi cua 1 mDl
6 *Cp; co Hi*3 la Cpl tioaç Cpu
7 vcpl : Sìlìaat đunự cua Chat lo n ;
« *L'pv;rtfiiet dưng cua pha hoi
9 sw itch di
10 case ('ioncr qua la n h ')
11 q - 1 + C p*(c - t£ ) /la f ld a ( x ,h v ) ;
12 Sct:o3U5t <90 941
13 case ( 'dien. 301 ' )
14 ct ■ 1;
15 ỡ*Sẽ ('m ot phan h ò i')
16 q; = ŨL;
1? case ( 'd ie u ngunạ t u ')
13 C£ » o^-
19 case t'n h ie c qua s o i ' 5
20 q - C p*(ĩ - tỄ ) /la a d a (K jh v );
21 í u n í i i e t t a i diem. nọutig tu
22 ù th e tu is e
23 d is p ( 1b ie n diet; k 1er. chi duoc chap nhan i tro n g 4 tru o n g hop t r e n 1)
24 end
25 q_3lope = (Ị/(q -l) ;
2« *he so goc
27
lamtìí.m Srtloirií.m material.m q Hne.m

q line Ln 1 Oal 1
Hình 7.6
Ví du : dk = ‘long qua lanh’ (nguyên liệu nạp vào là lỏng và ở nhiệt độ thấp),
t =368 K , t f =293 K .
[ q , q_slope] = q_line(Ểlong qua lành’, 368 ,293) = [1.3715 3.6917]
Nghĩa là : q =1.3715 và q_slope = 3.6917
Muc đỉch : tính giá trị của q trong điều kiện chưng cất và hệ số góc của đường q.
5 .3.9 Hàm ‘vpequiLm’
C ode:
119
'vr:\M ATlAB\hieuVdiJ1¡1IVvpequií.m
1Fte Ed*Tex; W refcvi H efc.
□ y m m M.
functiony-vptquil(v>
ẹlobai zK 5UB=0-
n=length.(K);
ío t i= l: n
tu . <KU)-l>*Stt>;
m au»l+U(K(i)-l)*v;
sun » 3 B+ tu/m au;
end
y=sua;
la m d a m antoine m Tialenal m q _ lin e /n vpéQuilm
vpeauil Ln 10 c 0 1 1Ệ
Hình 5.7
Cú pháp : y = ypequil(v)
%y: là hiệu của tổng cảc thành phần pha hơi với tổng các thành
phần pha lỏng cân bằng với nó
%v : tỷ số của lượng pha hơi với lượng pha lỏng
Vỉ du : xem bài tập 2 trong mục 3.1.2
Muc đich : giải phương trình vpequil(v) = 0 để tìm giá trị V. Sở dĩ có điều này là
do sự kết hợp của phương trình 5.7 và 5.10 . tò giá trị của V đó ta có thể xác định được
thành phần các cấu tử trong pha lông cũng như trong phã hơi cân bằng với nó.
5.4.1 Hàm ‘ActEn.m’
C ode:
function E = ActEn(T)
global k
R = 8.314;
%cong thuc tỉnh nang luong ho at hoa
E = log(k(2)/k(l ))*R*(T(2)*T( 1)/(T(2)-T(1)));
Củ pháp : Ea = ActEn(T)
Ví du : cho k = [ 1.6 2.8] X104, Ea = ActEn([540 550]) = 311308 J/mol
Minh họa trong hình 7.8 .
120
Muc đích : tính năng lượng hoạt hóa của phản ứng theo phương trình Arrheneous.
5.4.2 Hàm ‘Arrheneous.m’
C ode:
function y = Arrhenous(T)
global E kO
R = 8.314;
y = kO*exp(-E/(R*T));
Cứ pháp : k = Arrheneous(T)
Muc đích : tính hằng số tốc độ phẳn úng ở nhiệt độ nhất đỉnh khỉ biết hệ số nhân
kO, năng lượng hoạt hóa Ea của phản ủng.
He tc»t Vietv Web Window Help
D GÍ CỄ % ọ Currenỉ Drertory EYMTLABHeưtimAcs

□
ma
E:\HAT V Q Ẽ Cf M » g lo b a l k
» k • [ 1 . ỉ , 9*4] ;
Al l F i le s F ile » Ea = À cĩE n t[5 4 í 550 ]>
^ A C T In . l H-Í1
Eo -
[ 3 A iihensoư s. a H -fij
P53RJí.n H -íi' 1.39l8e+ 00S
Ị ^ b r l.n H“f i
Ìb L 2 .. H - tí »
[ 3 b ĩ 3.0 Hfj
[3 iM .l H -fl!
O lH F x .il n - fi!
p$ Sncp3.jpg JPC ỉ
3 to .a H -íỉ
1< ||«3___________ >
Current DlrectủỊV
-Ạ s ta rt'I
Hình 5.8
5.4.3 Hàm ‘BRk.m’
C ode: như ữong nội dung hình 5.9.
% x a : độ chuyển hóa của phản ứng
% t: thời gian mà phản úng đạt được đến độ chuyển hóa xa
%n : bậc của phản ứng ( ở đây tôi mới chỉ thiết lập cho phản ứng bậc 1
và bậc 2)
% const: tỷ lệ của hệ số phản ứng (A+2B -> C)thì const = —,Giả sử cho
Tùy thuộc vào phản ứng mà chỉnh const cho phù hợp
Ky thuat tinh toin trong cong nghp h6a hpc
%quad8: hàm tính tích phân

Củ pháp : k = BRk(xa)
E:\MATLAB\hieu\kirieticsWlkm (3 Ẽ ®
File Edit Text Window Help
□ Ẽ? % m V d M (2)
1 fu n c tio n y = BRkJxaJ
2 g lo b a l t n co n st
3 i f n ==1
4 S = quads(@ ylr 0#x a , le - 4 J / t ;
5 elsexE n ==2
6
1 y = wadS (@y2,0xx a , l- 5 - 4 ) / t ;
0 end
9 fu n c tio n y = y l(K )
10 Y = l. / ( l- x ) ;
11 fu n c tio n y = y2<x)
12 g lo b a l co n st
13 Y = l . / t (JL-jc) . * ( 1 - con st. * x ) ) ;
ActEn.m A rrh en e au s.m BRk.m

BRk Ln 1 con
Hình 5.9
Muc ãich : tính hằng sổ cân bằng của phản ứng trong bình phản ứng gián đoạn
theo phương trình 4.21
Ví du : cho phản ứng là bậc 1 , thời gian lưu là 0.5 h , và độ chuyển hóa là 0.75 .
Ta có ; k = BRk(0,75) = 2.7726 . Minh họa như hình 7.10.
• i MAILMi E PP«
I f is rar viaM w* Vitrine »fell
1 ũ ĩ£
lü Q
E: \1UT - Q (t) io ợl^lol t n const
» K • 0.3:
A il m e i n ie » n - ii
S l c u * .* B-Í1. » k - BWIÖ.7S)
n-ta Bacniwf: «ỊUAD0 18 obaoLc-br. qtUDL I» lUi 1 -c.w r.V J iqtl«ciuDti
¡Type "coznlng »fC HJßTfcS:«ĩu<iđ8:ữb50lcccfũncti0Ji‘‘ to suppress this
Q srk. uv ẰS? ; > ¿ti C:\IWTUD8^5VMỞlbởxV^ầtlẽb^f t * f t « »e *|A* 35
®BRk.B H-Ci. la C: ’ilULTUUB'ihicuUclatlicSV.BRkji at line 4
B1 n-fi.
k-
[jjM .si B-Ỉ1.
S b t 3 .t n-Ểí 2.TĨ2Ễ
■¿¿mw.t. n-n
1%IHR«. 1 E 21.
[;]flnapz.»g JJG ;
jjp r/i. a H-Í1
las >
1 Current t>ir?Uwy
Hình 5.10
122
5.4.4 Hàm ‘IMRx.m’

C ode:
function y = EMRx(xa)
global caO cbO n to k
if n = = 1
y = xa/(k*(l-xa))~to;
eiseif n = =2
y = xa/(k*cbO*(l -xa)*(l -caO*xa/cbO))-to;
end
Muc đích : giải phương trình IMRx(xa) = 0 để tìm độ chuyển hóa của sản phẩm
khỉ ra khỏi thiết bị phản ứng liên tụ.c ừong thời gian lưu biểu kiến cho trước.
Vỉ du : xem bài tập 3 trong mục 4.2,3
5.4.5 Hàm ‘to.m’
Code : minh họa trong hình 7.11
-.-1 E:\MATLABVhieu\kinetics\to.m B E D
I File Edit Text Window Help
I □ fl u % ü m J4 / , d
1 fw ncfcion y = to(d k jX aJ
2 g lo b al k n ep const
3 if n ==1
4 if dk == ' SIR'
5 y = q u a . d 8 ( @ ¥ l f Qlf:x:a1, l e - 4 j , / k ;
5 e l s e a £ dfc == ' IHR!
7 ¥ = £xa) / h ;
8 end
9 e lse if n ==2
10 if cte == 1PK R 1
11 Y = q u a d S (@y2 f O j - x a j l e - 4 ! y k ;
12 e lse if dk = = 'IM R '
13 ¥ 4 (x a)/£ ;
14 end
15 end
16 fu n ctio n y = Y-Ux)
17 g lo b al ep
15 5? = ( l+ e p .'* x ) . / ( 1 - x ) ;
19 fu n ctio iii y = y2(x)
20 g lo b al c o n s t ep
21 Y = (1+ ep. * x ) . J [ < l-x j • * (1- c c n s t.'* « ) ) ;
22 fu n ctio n y = y3(x )
23 g lo b al ep
24 2 = k . * ( 1 + e p . ’' x } . / ( 1 - x J ;
25 fu n ctio n y = T ^ lx )
25 g lo b al c o n s t ep
27 5 = x . * (1-hep.i x ) . / ( ( 1 - x ) .* (L- c o n st.* x ));
to Ln 1 Col 1
Hình 5.11
123
%k : hằng số tổc độ phản ứng

%n : bậc phản ứng
%ep : hệ số giãn nở thể tích (£•)
% const: giống như của hàm BRk.m
Cú vhảp : t = to(dk, xa)
% t : thời gian lưu biểu kiến của thiết bị làm việc , có thể là thiết bi
khuấy lỷ tưởng liên tục hoặc có Ihể là thiết bị ống dòng lý tưởng.
%dk : thiết bị đang khảo sát , ống dòng lý tưởng (PFR) hay thiết bị
khuấy lý tưởng liên tạc (IMR)
%xa : độ chuyển hóa của sản phẩm ở đầu ra của thiết bị
Muc đich : tính thời giãn lưu biểu kiến của các thiết bị IMR hãy thiết bị PFR.
Vỉ du : cho bậc của phân ứng là bậc 1 , hằng số tốc độ phản ứng k = 5.375 h “1 ,
độ chuyển hóa là 0.75 , hệ số giãn nở thể tích e = 2A. Thì ta có :
t l = to(‘PFR’, 0.75) = 0.542 h
t2 = toCIME*.’ , 0.75) = 1.5628 h
Minh họa trong hình 5.12.

ĐE®
File Ec* V w Wefc Wirriow Help
ũ |¿ c ft tí- ? Cwftwil H ríc ta y . it-'MAlLttirtTieurJiirwKS
r]Q
SH
E: VEiT V Q 0 ữ ít glób&L k n ep consc
3Ẽ = 5- 375*"
A i: F ile s f ile
U à c tE n .u H -fỉ cp - 2,4;
Acrhenéuu?.» K-M I» t l - t e ( 'í ,n«, ,C,'?5)
3 BP*■a K - t l|
R-11
A -É ll ũ. í 420
^ b t 3.m M -fi
tỉ ■ to riH R ,C.75J
^ b c 4 .n K -£ i
K t i.
[;J 5nap’2 ripc oro :
3 to .n M -fi
< Í 1 ___________ > -H
C u rren! □ itỂ tlu iy
•ặiaart
Hình 5.12
124
Tài liệu tham khảo
1) Dilip K. Das , P. E. , Dr. Rajaram K. Prabhudesai, P. E. , EIT Chemical

Review, 2s edn . Engineering Press : USA, 1998
2) Bruce A. Finlayson , Introduction to chemical engineering computing. John
Wilay & Sons , Inc., 2006
3) Ela P, ekalska, Maijolein van der Glas , Introduction to Matlab. Pattern
Recognition Group, Faculty o f Applied Sciences Delft University o f Technology ,
January 2002
4)
http://www.mathworks.com/matlabcentral/link exchange/MATLAB/Academic Curri
cula/Engineering/Chemical Engineering/
5) http://www.che.ncsu.edu/current students/matlab che.pdf
6) http://www.owlnet.rice.edu/~chbe403/bdistmccabe.html
125

Kỹ thuật tính toán trong công nghệ hóa học PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Kỹ thuật tính toán trong công nghệ hóa học PDF

Uploaded by

Copyright:

Available Formats

Kỹ thuật tính toán ừong công nghệ hóa học

Chương 8: Văn bân.................................................................................................................43

Lời nói đầu

Chương I:Bắt đầu với Matỉab

1.1 Nhập dữ liệu qua dòng lệnh

1.2 Sử dụng help trực tuyến

1.3 Đường dẫn

1.4 Lưu và tải dữ liệu

Chương 2: Các cấu trúc cơ bản và biến

2.1 Tính toán với Matlab

2.2 Giới thiệu dạng dữ liệu

Chương 3: Tính toán vói matrận và vec tơ

3.1.1 Kỉ tự hai chẩm và ừích các phần tử của vec tơ

» r(3:6) % lẩy ra các phần từ của r cổ vị tri từ 3 đến 6

>> r(5:-l:2) % chuyện gì sẽ xảy ra với trường hợp này?

3.1.2 Vec tơ cột và sự chuyầĩ vị

3.1.3 Nhân, chia và mũ của vec tơ

Một cách khác để tính tích vô hướng là sử dụng dấu

??? Error using ==> /

3.2.1 Các ma trận đặc biệt

» A(4,3) % this is not possible: A is a 3-by-3 matrixl

3.2.3 Các phép toán với ma trận

> > B-2 % ừ ừ các phần tử của B cho 2

4.1 Đồ thị đom giản

Mâu của dường vẽ Ký tụ hiên thị Nét vẽ

Hình J_i ĐỒ thị hàm số f i x ) = X3

Hình 3: vẽ đồ thị nhiều hàm trên cùng một hình

> > z=(x. A2). *y-2 *y;

Hình 5: Vẽ mặt f ( x ty) = x 2y - 2y ữong 3 chiều

/ (x, y) = . ^ trên lưới vuông [-1,1] X[-1,1] • Thực hiện lệnh

Hình 6a: Vẽ măt fipc, y) = . xy

Hình 6b: Vẽ mặt ừong 3 chiều tránh khuyết

Chương 5: Các ỉuồng điều khiển

5.1 Các toán tử logic

Lệnh Kêt quả

A=(b~=c) A=1 nêu b khác c

A=(b&c) A=1 nêu b =true và c=true

5.2 Lệnh find

>> plot (x,y, ’) ; hold on

Chúng ta cũng có thể xác định vị trí theo hàng cột.

» [I.J] =fin d (A>= 2.5)

*)If... else ...end

*) I f ... elseif... else ... end.

_____________ Syntax __________ ______ Exaunple ____

_______________ S y n t a x ____________ _______________ E xa m p le ____________

5.4 Cấu trúc lặp

Chỉ Số của phép lặp có thể là vec tơ:

for i-t:2:n for i-n:-l:3 for 1-0:0.6:4 for X - [26 9 81]

*) Vòng lặp while

Chương 6: Phương pháp sổ

6.1 Đường cong hồi quy

6.2 Phép nội suy

6.3 Giá trị của hàm sỗ

6.3.1 Hàm 'inline’

6.4 Phép tính tích phân và tích phân

6.5 Tính toán sổ học và cấu trúc điều khiển

while (x > le-6)

Sẽ cho ta kết quả gần đúng của hàm / = — .

Chương 7: Viết chương trình trong Matlab

7.1 Kịch bản và hàm m-file

7.2 Hàm m-file

function [avr.sd] = stat(x)

7.2.1 Những biến đặc biệt của hàm số

> > z=(x. A2). y-2 y;

_____ Syntax Exaunple

_______ S y n t a x _ E xa m p le ______

Nu, = 0.41 5745060.687°33 = 65.2 => 03= ^ ^ = ^ ^ ^ ^ - = 98.5 W/m2K