Bước hoàn thiện trong quy trình thực nghiệm xác định hằng số động học là tính hệ số sinh khối Y và hằng số phân hủy nội sinh kp từ phương trình PT 6-34: 1/θc = Y·(So-S)/(τ·X) – kp. Phương trình này là đường thẳng với vsu/X làm biến số, cho phép xác định Y (độ dốc dương) và kp (trị tuyệt đối của tung độ gốc âm) bằng hồi quy tuyến tính. Bài viết trình bày đầy đủ phương pháp bình phương cực tiểu, ví dụ số hoàn chỉnh với τ = 6 giờ và 5 mức θc, so sánh tuyến tính hóa đơn giản với tối ưu hóa số phi tuyến, tiêu chí kiểm tra chất lượng bộ hằng số động học, và hướng dẫn ứng dụng thực tiễn cho các loại nước thải đặc trưng tại Việt Nam.
1. Nền tảng toán học của PT 6-34
1.1 Xuất xứ và cấu trúc
Từ phương trình thăng bằng sinh khối ở trạng thái ổn định (PT 6-29): 1/θc = Y·|vsu|/X – kp, thay |vsu| = (So-S)/τ:
1/θc = Y · (So-S)/(τ·X) – kp — (PT 6-34)
Đặt y = 1/θc (trục tung) và x = (So-S)/(τ·X) (trục hoành), đây là đường thẳng y = Y·x – kp với độ dốc = Y (dương) và tung độ gốc = -kp (âm).
1.2 Ý nghĩa vật lý của độ dốc Y
Độ dốc Y = Δ(1/θc) / Δ(vsu/X) biểu thị: khi tốc độ tiêu thụ cơ chất riêng tăng thêm 1 mg COD/(mg VSS·ngày), tốc độ cần thải bùn tăng thêm Y ngày⁻¹ để duy trì X ổn định. Điều này phản ánh trực tiếp: Y mg sinh khối mới được tạo ra từ 1 mg COD tiêu thụ, và lượng sinh khối mới đó cần được thải ra.
1.3 Ý nghĩa vật lý của tung độ gốc -kp
Khi vsu/X = 0 (không có cơ chất), PT 6-34 cho 1/θc = -kp < 0, tức θc âm — không có nghĩa vật lý. Thực tế, điều này chỉ là phép ngoại suy đường thẳng ra ngoài vùng dữ liệu. Ý nghĩa đúng là: điểm cân bằng tối thiểu là vsu/X = kp/Y — tốc độ tiêu thụ cơ chất tối thiểu để bù đắp phân hủy nội sinh. Nếu vsu/X < kp/Y, sinh khối trong hệ sẽ suy giảm dần.
Với Y = 0,5 và kp = 0,06 ngày⁻¹: vsu/X_min = 0,06/0,5 = 0,12 mg COD/(mg VSS·ngày). Nếu tốc độ tiêu thụ cơ chất riêng rơi xuống dưới mức này, hệ mất dần sinh khối — thực tế xảy ra khi bể bị sự cố đói cơ chất (thiếu nước thải đầu vào).
2. Phương pháp bình phương cực tiểu — Tính Y và kp
2.1 Công thức tính toán
Với n cặp số liệu (xi, yi), phương trình hồi quy y = Y·x + b (b = -kp):
Y (độ dốc) = [n·Σxi·yi – Σxi·Σyi] / [n·Σxi² – (Σxi)²]
kp = -b = -(ȳ – Y·x̄) = Y·x̄ – ȳ
Trong đó x̄ = Σxi/n, ȳ = Σyi/n.
2.2 Ví dụ tính toán — Nước thải sinh hoạt đô thị, Hà Nội
Điều kiện thực nghiệm: V = 10 L, Q = 40 L/ngày → τ = V/Q = 0,25 ngày (6 giờ), So = 250 mg COD/L, T = 20°C.
Dữ liệu thực nghiệm tại 5 mức θc ổn định:
| θc (ngày) | S (mg COD/L) | X (mg VSS/L) |
|---|---|---|
| 5 | 25 | 1.500 |
| 8 | 15 | 2.200 |
| 12 | 10 | 2.900 |
| 16 | 8 | 3.400 |
| 20 | 7 | 3.800 |
Tính các cột hỗ trợ:
| θc | x = (So-S)/(τ·X) | y = 1/θc | x² | x·y |
|---|---|---|---|---|
| 5 | (250-25)/(0,25×1500) = 0,600 | 0,200 | 0,3600 | 0,12000 |
| 8 | (250-15)/(0,25×2200) = 0,427 | 0,125 | 0,1823 | 0,05338 |
| 12 | (250-10)/(0,25×2900) = 0,331 | 0,083 | 0,1096 | 0,02747 |
| 16 | (250-8)/(0,25×3400) = 0,285 | 0,063 | 0,0812 | 0,01796 |
| 20 | (250-7)/(0,25×3800) = 0,256 | 0,050 | 0,0655 | 0,01280 |
| Σ | 1,899 | 0,521 | 0,7986 | 0,2316 |
Áp dụng công thức bình phương cực tiểu (n = 5):
Y = (5 × 0,2316 – 1,899 × 0,521) / (5 × 0,7986 – 1,899²) = (1,158 – 0,9894) / (3,993 – 3,6062) = 0,1686 / 0,3868 = 0,436 g VSS/g COD
x̄ = 1,899/5 = 0,3798; ȳ = 0,521/5 = 0,1042
kp = Y·x̄ – ȳ = 0,436 × 0,3798 – 0,1042 = 0,1656 – 0,1042 = 0,0614 ngày⁻¹
Kết quả: Y = 0,44 g VSS/g COD, kp = 0,061 ngày⁻¹
Kiểm tra hợp lệ: Y = 0,44 nằm trong khoảng 0,4-0,6 cho nước thải sinh hoạt (hợp lý). kp = 0,061 ngày⁻¹ nằm trong khoảng 0,05-0,10 ngày⁻¹ (hợp lý).
3. Kiểm tra chất lượng bộ số liệu
3.1 Hệ số tương quan R²
Tính R² để đánh giá mức độ phù hợp của đường thẳng hồi quy:
Với đường hồi quy ŷi = 0,436·xi – 0,0614:
| θc | yi (thực) | ŷi (tính) | (yi-ŷi)² | (yi-ȳ)² |
|---|---|---|---|---|
| 5 | 0,200 | 0,436×0,600-0,061=0,201 | 0,000001 | 0,009177 |
| 8 | 0,125 | 0,436×0,427-0,061=0,125 | 0,000000 | 0,000432 |
| 12 | 0,083 | 0,436×0,331-0,061=0,083 | 0,000000 | 0,000449 |
| 16 | 0,063 | 0,436×0,285-0,061=0,063 | 0,000000 | 0,001697 |
| 20 | 0,050 | 0,436×0,256-0,061=0,051 | 0,000001 | 0,002937 |
| Σ | SSres=0,000002 | SStot=0,014692 |
R² = 1 – SSres/SStot = 1 – 0,000002/0,014692 ≈ 0,9999
R² xấp xỉ 1,0 cho thấy số liệu ví dụ này được tạo ra rất “sạch” (ít sai số). Trong thực tế, R² = 0,95-0,99 là kết quả tốt, R² = 0,90-0,95 là chấp nhận được.
3.2 Kiểm tra giá trị kp/Y — Nồng độ cơ chất tối thiểu
Nồng độ cơ chất tối thiểu để hệ tự duy trì được tính từ PT 6-12 khi μ’ = 0:
μ’min = 0 = k·S_min/(Ks+S_min) – kp
S_min = Ks·kp / (k – kp)
Với k = μm/Y (hằng số tốc độ tiêu thụ cơ chất), Ks, Y, kp từ thực nghiệm.
Ví dụ: nếu k = 0,30 mg COD/(mg VSS·ngày), Ks = 40 mg COD/L, kp = 0,061 ngày⁻¹:
S_min = 40 × 0,061 / (0,30 – 0,061) = 2,44 / 0,239 = 10,2 mg COD/L
Nghĩa là hệ bùn hoạt tính với các hằng số này chỉ có thể xử lý COD xuống mức tối thiểu khoảng 10 mg/L trong điều kiện lý tưởng. Đây là thông tin quan trọng cho kỹ sư đánh giá giới hạn xử lý của nước thải cụ thể.
4. So sánh phương pháp tuyến tính hóa đơn giản và tối ưu hóa bình phương cực tiểu phi tuyến
4.1 Phương pháp tuyến tính hóa đơn giản (PT 6-33 và 6-34)
GS. Lê Văn Cát ghi nhận: “Phương pháp xác định các hằng số tốc độ theo phương pháp tuyến tính hóa đơn giản dễ thực hiện nhưng cho độ chính xác về mặt tính toán không cao.”
Ưu điểm:
- Tính toán bằng tay hoặc Excel đơn giản
- Không cần phần mềm chuyên dụng
- Dễ kiểm tra và giải thích kết quả
- Phù hợp cho thực hành kỹ thuật thông thường
Nhược điểm:
- Biến đổi nghịch đảo khuếch đại sai số tại S nhỏ (điểm S thấp, 1/S cao)
- R² của đường thẳng nghịch đảo không phản ánh trực tiếp chất lượng khớp mô hình Monod gốc
- Khi sai số đo lớn (thường trong thực hành), kết quả k, Ks có thể dao động nhiều
4.2 Phương pháp tối ưu hóa bình phương cực tiểu phi tuyến
“Cũng từ các số liệu trên có thể tính các thông số động học thông qua phương pháp tối ưu dựa trên cơ sở nguyên tắc bình phương cực tiểu sai số của số liệu thực nghiệm và số liệu tính toán.” — GS. Lê Văn Cát.
Phương pháp này tìm k, Ks, Y, kp sao cho tổng bình phương sai số giữa S thực nghiệm và S tính từ mô hình Monod là nhỏ nhất, không cần biến đổi phi tuyến.
Hàm mục tiêu tối thiểu hóa:
F(k, Ks, Y, kp) = Σᵢ [S_i,thực – S_i,tính(k, Ks, Y, kp)]²
S_i,tính được tính bằng cách giải hệ phương trình trạng thái ổn định cho bể CSTR với hồi lưu.
Ưu điểm:
- Không có méo sai số do biến đổi nghịch đảo
- Tối ưu hóa trực tiếp trên không gian mô hình gốc
- Cho kết quả chính xác hơn đặc biệt khi số liệu có nhiều sai số
- Có thể tối ưu hóa đồng thời tất cả 4 hằng số cùng lúc
Nhược điểm:
- Cần phần mềm tối ưu hóa phi tuyến (MATLAB, Python scipy.optimize, R package)
- Có thể hội tụ về cực trị địa phương nếu điểm xuất phát không tốt
- Phức tạp hơn và khó giải thích trực quan
4.3 Khi nào dùng phương pháp nào?
Trong thực hành tại Việt Nam, khuyến nghị:
Dùng tuyến tính hóa (PT 6-33, 6-34) khi:
- Số điểm dữ liệu ≥ 6 và phân bố tốt
- Sai số đo nhỏ (CV < 10%)
- Cần kết quả nhanh cho báo cáo nghiên cứu sơ bộ
- Không có phần mềm tối ưu hóa phi tuyến
Dùng tối ưu hóa phi tuyến khi:
- Thiết kế chính thức công trình xử lý nước thải quy mô lớn (> 10.000 m³/ngày)
- Số liệu có sai số cao (môi trường thí nghiệm không lý tưởng)
- Muốn ước tính khoảng tin cậy (confidence interval) của hằng số
- Sử dụng mô hình ASM phức tạp hơn (ASM1, ASM2, ASM3)
Thực hành tốt nhất: chạy cả hai phương pháp và so sánh kết quả. Nếu kết quả sai lệch < 15%, có thể tin tưởng dùng kết quả tuyến tính hóa. Nếu sai lệch > 20%, cần dùng tối ưu hóa phi tuyến.
5. Yêu cầu chất lượng thực nghiệm và tiêu chí hợp lệ
5.1 Yêu cầu về số lượng điểm dữ liệu
GS. Lê Văn Cát quy định rõ: “thí nghiệm cần được tiến hành với S biến động trong khoảng rộng, tức là với giá trị θc rất khác nhau (không ít hơn 5 giá trị) với trị số từ 1 – 10 ngày.”
Tuy nhiên, trong thực tế:
- 5 điểm là tối thiểu — đủ để vẽ đường thẳng có ý nghĩa thống kê
- 7-8 điểm là khuyến nghị — giảm ảnh hưởng của điểm bất thường
- 10+ điểm là lý tưởng — cho phép phát hiện quan hệ phi tuyến nếu có
Phân bổ θc nên đều nhau trên thang logarithm: ví dụ θc = 2, 3, 5, 8, 12, 18 ngày (thay vì 2, 4, 6, 8, 10 ngày). Điều này đảm bảo S = S(θc) phân bổ đều hơn, tránh tập trung quá nhiều điểm ở vùng S cao.
5.2 Tiêu chí chấp nhận bộ hằng số động học
Sau khi tính được k, Ks, Y, kp, kiểm tra theo bảng tiêu chí sau:
| Hằng số | Khoảng hợp lệ (nước thải sinh hoạt, 20°C) | Lý do từ chối nếu ngoài khoảng |
|---|---|---|
| k | 0,10-0,40 g COD/(g VSS·ngày) | Sai số đo X hoặc S |
| Ks | 20-150 mg COD/L | Số điểm S thấp không đủ đại diện |
| Y | 0,35-0,65 g VSS/g COD | Sai số đo X (VSS) |
| kp | 0,04-0,12 ngày⁻¹ | Sai số tính θc hoặc không đạt trạng thái ổn định |
5.3 Kiểm tra nhất quán nội tại
Từ bốn hằng số, kiểm tra tính nhất quán nội tại:
- μm = k·Y nên nằm trong 0,05-0,20 ngày⁻¹ cho vi khuẩn dị dưỡng hiếu khí ưa nhiệt vừa.
- S_min = Ks·kp/(k-kp) nên dương và < 20 mg COD/L (hệ bùn hoạt tính có thể xử lý xuống nồng độ rất thấp).
- SRT tối thiểu = 1/(k·Y·So/(Ks+So) – kp) với So điển hình nên < 2 ngày (có thể khởi động hệ nhanh).
- Dự đoán S tại từng θc bằng mô hình Monod và so sánh với S thực nghiệm — sai số < 15% cho mọi điểm.
6. Ứng dụng bộ hằng số động học trong thiết kế công trình
6.1 Từ thực nghiệm đến thiết kế — Chuỗi tính toán
Sau khi có đủ bốn hằng số, kỹ sư tiến hành thiết kế theo trình tự:
Bước 1 — Hiệu chỉnh nhiệt độ:
- k_T = k_20 × (1,09)^(T-20) (PT 6-24)
- Ks_T = Ks_20 × (1,09)^(T-20)
- kp_T = kp_20 × (1,04)^(T-20) (PT 6-23)
- Y: không hiệu chỉnh đáng kể
Bước 2 — Tính SRT thiết kế: Chọn SRT sao cho nồng độ S đầu ra đạt tiêu chuẩn, với hệ số an toàn Fs = 1,5-2,5.
Bước 3 — Tính thể tích bể (V): V = Q × SRT × (Y × (So-S)/X – kp) … rút từ phương trình cân bằng sinh khối
Bước 4 — Tính lượng bùn sản sinh: Px = Y’ × Q × (So-S) / 1000 (kg VSS/ngày)
Y’ = Y / (1 + kp × SRT) (hệ số sinh khối thực tế)
Bước 5 — Tính nhu cầu oxy: O2 = Q × (So-S) × (1 – 1,42 × Y’) / 1000 (kg O2/ngày)
6.2 Ví dụ ứng dụng cho nước thải khu công nghiệp thực phẩm
Một khu công nghiệp chế biến thực phẩm tại TP.HCM có:
- Q = 500 m³/ngày
- So = 800 mg COD/L
- T = 30°C (miền Nam)
- Tiêu chuẩn xả thải: S ≤ 50 mg COD/L
Hằng số động học từ thực nghiệm pilot (ở 20°C):
- k_20 = 0,28 mg COD/(mg VSS·ngày), Ks_20 = 100 mg COD/L
- Y = 0,55 g VSS/g COD, kp_20 = 0,07 ngày⁻¹
Hiệu chỉnh nhiệt độ về 30°C:
- k_30 = 0,28 × (1,09)^10 = 0,28 × 2,367 = 0,663 mg COD/(mg VSS·ngày)
- Ks_30 = 100 × 2,367 = 237 mg COD/L
- kp_30 = 0,07 × (1,04)^10 = 0,07 × 1,480 = 0,104 ngày⁻¹
Chọn SRT = 10 ngày. Tính X từ phương trình cân bằng:
X = Y × (So-S) × SRT / (τ × (1 + kp × SRT))
Đặt τ = 0,25 ngày (6 giờ), S ≈ 50 mg COD/L (mục tiêu):
X = 0,55 × (800-50) × 10 / (0,25 × (1 + 0,104 × 10)) = 0,55 × 750 × 10 / (0,25 × 2,04) = 4125 / 0,51 = 8.088 mg VSS/L — quá cao, cần tăng τ
Thử τ = 1 ngày (24 giờ): X = 4125 / 2,04 = 2.022 mg VSS/L (hợp lý, nằm trong 1.500-4.000 mg/L)
Thể tích bể: V = Q × τ = 500 × 1 = 500 m³
Lượng bùn sản sinh: Y’ = 0,55 / (1 + 0,104 × 10) = 0,55 / 2,04 = 0,270 g VSS/g COD Px = 0,270 × 500 × (800-50) / 1000 = 0,270 × 500 × 0,750 = 101 kg VSS/ngày
Nhu cầu oxy: O2 = 500 × (800-50) × (1 – 1,42 × 0,270) / 1000 = 500 × 0,750 × 0,617 = 231 kg O2/ngày
Kết quả thiết kế: bể CSTR V = 500 m³, X ≈ 2.000 mg/L, Px = 101 kg/ngày, O2 = 231 kg/ngày — đây là đầu vào trực tiếp cho thiết kế các hạng mục máy sục khí, bể nén bùn, máy ép bùn.
7. Sai số phổ biến trong thực nghiệm Việt Nam và cách khắc phục
7.1 Lấy mẫu khi chưa đạt trạng thái ổn định
Đây là lỗi phổ biến nhất, xảy ra khi kỹ sư muốn rút ngắn thời gian thực nghiệm. Biểu hiện: điểm dữ liệu nằm lộn xộn trên đồ thị, R² thấp < 0,85.
Khắc phục: vận hành ít nhất 3×θc ngày liên tục ở mỗi mức trước khi lấy mẫu. Với θc = 10 ngày, cần 30 ngày vận hành liên tục. Đây là yêu cầu không thể rút ngắn.
7.2 Phân tích VSS không chuẩn
VSS (Volatile Suspended Solids) là thước đo nồng độ vi sinh. Sai số đo VSS lớn (> 15%) dẫn đến sai số Y và kp tương ứng. Nguyên nhân: nhiệt độ nung không đúng (550°C ± 50°C), thời gian nung không đủ (2 giờ), giấy lọc không tiêu chuẩn.
Khắc phục: hiệu chuẩn lò nung, dùng giấy lọc GF/C hoặc tương đương, thực hiện mẫu trắng (blank), lặp lại ít nhất 3 lần mỗi điểm.
7.3 Bình phản ứng không khuấy trộn đều
Nếu có vùng chết (dead zone) hoặc gradient nồng độ trong bình, điều kiện CSTR lý tưởng không đạt, phương trình thăng bằng chất sai lệch. Kiểm tra bằng cách đo nồng độ S tại nhiều điểm trong bình — nếu S dao động > 10%, cần tăng tốc độ khuấy.
8. Kết luận
PT 6-34 cung cấp phương pháp đơn giản, hiệu quả để xác định hai hằng số động học còn lại Y và kp từ cùng bộ số liệu thực nghiệm đã dùng để tính k và Ks (PT 6-33). Hai phương trình đường thẳng — PT 6-33 cho cặp (k, Ks) và PT 6-34 cho cặp (Y, kp) — tạo thành một quy trình hoàn chỉnh xác định đủ bốn hằng số động học cốt lõi của mô hình Monod mở rộng.
Phương pháp tuyến tính hóa đơn giản (bình phương cực tiểu trên đồ thị thẳng) dễ thực hiện bằng Excel và đủ chính xác cho hầu hết ứng dụng kỹ thuật, trong khi phương pháp tối ưu hóa phi tuyến (số) cho kết quả chính xác hơn cho công trình quy mô lớn. Kết hợp cả hai phương pháp và kiểm tra chéo kết quả là thực hành tốt nhất để đảm bảo bộ hằng số động học đáng tin cậy — nền tảng vững chắc cho thiết kế công trình xử lý nước thải đạt tiêu chuẩn tại Việt Nam.