Nghiệm PT 7-26 Cho Bể Thứ 2 Và Công Thức Tổng Quát PT 7-27 Cho Bể Thứ i: Phân Bố Nồng Độ Dạng Erlang Trong Chuỗi CSTR Nối Tiếp
Giới Thiệu: Vẻ Đẹp Toán Học Của Phân Bố Nồng Độ Trong Chuỗi Phản Ứng
Khi giải phương trình vi phân PT 7-25 cho nồng độ C2 ở bể thứ 2 trong chuỗi n CSTR nối tiếp, GS. Lê Văn Cát nhận được kết quả:
C2 = C₀·n·H·e^(-n·H) (PT 7-26)
Và công thức tổng quát cho bể thứ i:
Ci = C₀/(i−1)! · (n·H)^(i-1) · e^(-n·H) (PT 7-27)
Hai kết quả này không chỉ là nghiệm toán học thuần túy — chúng có mối liên hệ sâu sắc với phân bố xác suất Erlang trong lý thuyết thống kê, và tiết lộ một đặc điểm hành vi thú vị: nồng độ tại mỗi bể trung gian không giảm đơn điệu từ C₀ về 0 ngay lập tức mà đạt cực đại tại một thời điểm xác định trước khi suy giảm — khác hoàn toàn với hành vi của bể đầu tiên.
Bài viết này phân tích chi tiết các nghiệm PT 7-26 và PT 7-27, ý nghĩa vật lý và toán học, đặc điểm cực đại, mối liên hệ với phân bố Erlang, và hàm ý thực tiễn cho thiết kế và vận hành hệ thống XLNT.
Phần 1: Dẫn Xuất Nghiệm PT 7-26 Cho Bể Thứ 2
1.1 Giải PT 7-25 Bằng Nhân Tử Tích Phân
Từ PT 7-25: dC2/dt + (n/τ)·C2 = (n/τ)·C₀·e^(-n·t/τ)
Nhân tử tích phân (integrating factor): μ(t) = e^((n/τ)·t)
Nhân cả hai vế với μ:
d[C2·e^(n·t/τ)]/dt = (n/τ)·C₀·e^(-n·t/τ)·e^(n·t/τ) = (n/τ)·C₀
Tích phân cả hai vế:
C2·e^(n·t/τ) = (n/τ)·C₀·t + A
Với điều kiện ban đầu t = 0, C2 = 0 (bể 2 ban đầu không có cơ chất — giả thiết chuẩn):
0 = 0 + A → A = 0
Vậy: C2 = (n/τ)·C₀·t·e^(-n·t/τ)
Viết theo H = t/τ:
C2 = C₀·(n·H)·e^(-n·H) = C₀·n·H·e^(-n·H) (PT 7-26)
1.2 Kiểm Tra PT 7-26 Tại Các Thời Điểm Quan Trọng
Tại H = 0: C2 = C₀·n·0·e⁰ = 0 ✓ (ban đầu bể 2 sạch)
Khi H → ∞: C2 = C₀·n·H·e^(-n·H) → C₀·n·∞·e^(-∞) → 0 ✓ (tiệm cận về 0 do e^(-n·H) giảm nhanh hơn n·H tăng)
Tại H = 1/n: C2 = C₀·n·(1/n)·e^(-n·1/n) = C₀·1·e⁻¹ = C₀/e ≈ 0,368·C₀
1.3 Đặc Điểm Nổi Bật: Cực Đại Của C2
Khác với C1 (giảm đơn điệu từ C₀), C2 tăng từ 0 đến cực đại rồi mới giảm về 0. Tìm cực đại bằng cách đạo hàm và cho bằng 0:
dC2/d(nH) = C₀·[e^(-nH) − nH·e^(-nH)] = C₀·e^(-nH)·[1 − nH] = 0
→ 1 − nH* = 0 → H = 1/n* (hoặc tương đương t* = τ/n = thời gian lưu một bể)
Giá trị cực đại: C2_max = C₀·n·(1/n)·e^(-n·1/n) = C₀·e⁻¹ ≈ 0,368·C₀
Lý giải vật lý: Lúc đầu (H nhỏ), nước thải từ bể 1 (giàu cơ chất C1 ≈ C₀) bắt đầu chảy vào bể 2 → C2 tăng. Sau thời điểm H = 1/n, nồng độ C1 đầu vào đã giảm đáng kể → lượng cơ chất nạp thêm vào bể 2 giảm dần trong khi phản ứng vẫn tiêu thụ → C2 giảm.
Phần 2: Công Thức Tổng Quát PT 7-27 Cho Bể Thứ i
2.1 Dạng Tổng Quát
Tương tự quá trình trên, bằng cách lặp lại giải ODE từng bước:
Ci = C₀/(i−1)! · (n·H)^(i-1) · e^(-n·H) (PT 7-27)
Trong đó:
- i: số thứ tự bể (i = 1, 2, 3, …, n)
- (i−1)!: giai thừa của (i−1)
- (n·H)^(i-1): lũy thừa (i−1) của n·H
- e^(-n·H): hàm mũ suy giảm chung
Kiểm tra với i = 1: C1 = C₀/0! · (n·H)^0 · e^(-n·H) = C₀·1·1·e^(-n·H) = C₀·e^(-n·H) ✓ (đúng với PT 7-24)
Kiểm tra với i = 2: C2 = C₀/1! · (n·H)^1 · e^(-n·H) = C₀·n·H·e^(-n·H) ✓ (đúng với PT 7-26)
Với i = 3: C3 = C₀/2! · (n·H)^2 · e^(-n·H) = (C₀/2)·(n·H)²·e^(-n·H)
2.2 Cực Đại Của Ci Tổng Quát
Đạo hàm PT 7-27 theo nH và cho bằng 0:
dCi/d(nH) = C₀/(i−1)! · [(i−1)·(nH)^(i-2)·e^(-nH) − (nH)^(i-1)·e^(-nH)] = 0
→ (i−1) − nH* = 0 → H = (i−1)/n*
Cực đại tại: t* = H*·τ = (i−1)·τ/n = (i−1) × (thời gian lưu một bể)
Bảng cực đại C_i theo i (với n = 3):
| Bể i | H* = (i-1)/n | t* (với τ = 12h) | Ci_max/C₀ |
|---|---|---|---|
| 1 | 0 | 0 | C₀ (giảm đơn điệu) |
| 2 | 1/3 = 0,333 | 4 giờ | C₀/e ≈ 0,368·C₀ |
| 3 | 2/3 = 0,667 | 8 giờ | C₀·4/(2·e²) ≈ 0,271·C₀ |
Ý nghĩa thực tiễn: Tại trạm XLNT khởi động với τ = 12 giờ, n = 3 bể:
- Bể 2 đạt nồng độ tối đa sau 4 giờ → đây là thời điểm tải trọng vi sinh ở bể 2 cao nhất
- Bể 3 đạt nồng độ tối đa sau 8 giờ → đỉnh tải sau bể 2 khoảng 4 giờ
Kỹ sư vận hành cần biết các thời điểm này để điều chỉnh lượng oxy cung cấp hoặc hóa chất trong giai đoạn khởi động.
Phần 3: Mối Liên Hệ Với Phân Bố Erlang — Toán Học Sâu Hơn
3.1 Phân Bố Erlang Là Gì?
Trong lý thuyết xác suất, phân bố Erlang bậc k mô tả tổng của k biến ngẫu nhiên độc lập, phân bố mũ với cùng tốc độ λ:
f(x; k, λ) = λ^k · x^(k-1) · e^(-λx) / (k-1)!
So sánh với PT 7-27 của GS. Lê Văn Cát:
Ci/C₀ = (n·H)^(i-1) · e^(-n·H) / (i-1)!
Đây chính xác là phân bố Erlang bậc i với tốc độ λ = n đánh giá tại x = H!
3.2 Ý Nghĩa Của Sự Tương Đồng
Tại sao phân bố nồng độ trong chuỗi CSTR lại là Erlang? Câu trả lời liên quan đến Phân Phối Thời Gian Lưu (RTD — Residence Time Distribution):
Trong chuỗi n CSTR nối tiếp, RTD của toàn hệ có dạng Erlang bậc n:
E(t) = (n/τ)·(n·t/τ)^(n-1)·e^(-n·t/τ) / (n-1)!
Đây chính xác là PT 7-27 với i = n và H = t/τ!
Ý nghĩa vật lý: Hàm RTD E(t) mô tả xác suất một phần tử chất lỏng lưu lại trong hệ đúng t giây. Khi n = 1 (CSTR đơn): Erlang bậc 1 = phân bố mũ E(t) = (1/τ)·e^(-t/τ). Khi n → ∞ (PFR): Erlang bậc ∞ → delta Dirac δ(t − τ) (tất cả lưu đúng τ giây).
3.3 Ứng Dụng RTD Trong Đánh Giá Thực Tế
Thí nghiệm dấu vết (tracer test) xác định RTD thực tế của công trình. Nếu RTD đo được phù hợp với Erlang bậc n, có thể kết luận:
- Hệ thống tương đương n bể CSTR nối tiếp về mặt thủy lực
- Số bể CSTR tương đương: n = τ² / σ²_RTD (trong đó σ² là phương sai của RTD)
- Nếu σ² < τ² → n > 1 → hệ gần PFR hơn CSTR đơn ✓
- Nếu σ² = τ² → n = 1 → hệ là CSTR đơn
Ví dụ thực tế: Đo RTD của bể aerotank hình chữ nhật tại Hà Nội, σ² = 0,25·τ²: n_tương đương = 1/0,25 = 4 bể CSTR tương đương
Điều này cho phép dùng PT 7-27 với n = 4 để mô phỏng hiệu suất xử lý thực tế của bể.
Phần 4: Đặc Điểm Đường Cong Ci Theo H — Phân Tích Định Lượng
4.1 Giá Trị Ci Tại Trạng Thái Ổn Định (H → ∞)
Từ PT 7-27, khi H → ∞: tất cả Ci → 0 (do e^(-n·H) → 0 nhanh hơn (nH)^(i-1) → ∞).
Điều này có nghĩa là trong hệ thống quá độ (bắt đầu từ bể sạch), nồng độ tại bể cuối cùng cuối cùng đạt 0 — mọi cơ chất đã bị tiêu thụ.
Nhưng đây là nghiệm của bài toán quá độ với điều kiện ban đầu C2 = 0. Khi bể được nạp liên tục, nghiệm thực tế tiệm cận về trạng thái ổn định C_i,∞ ≠ 0 (được tính theo PT 7-13 nhân chuỗi).
4.2 Bảng Giá Trị Ci/C₀ Theo H Cho n = 3
| H | C1/C₀ | C2/C₀ | C3/C₀ |
|---|---|---|---|
| 0 | 1,000 | 0 | 0 |
| 0,1 | 0,741 | 0,222 | 0,033 |
| 0,2 | 0,549 | 0,329 | 0,099 |
| 0,333 (H* C2) | 0,368 | 0,368 | 0,184 |
| 0,5 | 0,223 | 0,335 | 0,252 |
| 0,667 (H* C3) | 0,135 | 0,270 | 0,271 |
| 1,0 | 0,050 | 0,149 | 0,224 |
| 2,0 | 0,002 | 0,015 | 0,044 |
Nhận xét: C1 giảm đơn điệu; C2 tăng đến cực đại rồi giảm; C3 tăng muộn hơn, đạt đỉnh cao hơn (tương đối) và sau đó giảm. Đây là đặc trưng của “sóng nồng độ” truyền qua chuỗi bể.
4.3 Tốc Độ Lan Truyền Thông Tin Qua Chuỗi
Thời điểm bể i đạt nồng độ đáng kể (ví dụ Ci > 0,1·C₀):
Từ bảng, bể 3 đạt Ci > 0,1·C₀ khi nH > ~1,5, tức H > 0,5 (với n = 3).
Ý nghĩa vận hành: Sau khi khởi động hoặc sau sốc tải trọng, phải mất khoảng t = H*·τ = 0,5 × 12h = 6 giờ trước khi bể 3 (đầu ra) bị ảnh hưởng đáng kể. Đây là “thời gian trễ” (lag time) — cơ hội để kỹ sư vận hành điều chỉnh trước khi chất lượng đầu ra bị ảnh hưởng.
Phần 5: Kiểm Tra Bảo Toàn Khối Lượng
5.1 Tổng Cơ Chất Trong Toàn Hệ
Lượng cơ chất tổng trong n bể tại thời điểm H:
Tổng = (V/n)·∑[i=1 đến n] Ci
= (V/n)·C₀·∑[i=1 đến n] [(n·H)^(i-1)/(i-1)!·e^(-n·H)]
= (V/n)·C₀·e^(-n·H)·∑[i=1 đến n] (n·H)^(i-1)/(i-1)!
Khi n → ∞, ∑[i=1 đến ∞] (n·H)^(i-1)/(i-1)! = e^(n·H):
Tổng → (V/n)·C₀·e^(-n·H)·e^(n·H)·n/n = … (cần tính cẩn thận)
Trong giới hạn n → ∞ (PFR), tất cả cơ chất tập trung tại “wavefront” x = v·t = L·H và phân bố là hàm delta → bảo toàn khối lượng hoàn hảo.
5.2 Kiểm Tra Với n = 1 (CSTR Đơn)
Với n = 1, i = 1: C1 = C₀·e^(-H) = C₀·e^(-t/τ)
Đây chính xác là nghiệm PT 7-3 (phản ứng bậc nhất trong batch với t thay bằng H·τ = t). Nhất quán với các kết quả đã có.
Phần 6: Ứng Dụng PT 7-26, 7-27 Trong Phân Tích Sốc Tải Trọng
6.1 Phân Tích Sốc Tải Trọng — Khi C₀ Đột Ngột Tăng
Tình huống thực tế: Tại KCN chế biến thực phẩm Bình Dương, vào cuối ca sản xuất, COD đột ngột tăng gấp 3 lần trong 2 giờ trước khi trở lại bình thường. Hệ thống XLNT có n = 3 bể CSTR nối tiếp, τ = 18 giờ.
Dùng nguyên lý chồng chất (superposition): sốc tải trọng = tín hiệu bước nhảy + tín hiệu bước nhảy ngược sau 2 giờ.
Đáp ứng tại bể 3 (C3_sốc) = đáp ứng với C₀_mới − đáp ứng với C₀_mới (trễ 2 giờ)
Áp dụng PT 7-27 với i = 3, n = 3, τ = 18 giờ:
C3(t) = ΔC₀/(2!)·(3t/18)²·e^(-3t/18) – ΔC₀/(2!)·(3(t-2)/18)²·e^(-3(t-2)/18)·𝟙(t≥2)
Tính tại t = 8 giờ (khi H = 8/18 = 0,444):
- Số hạng 1: ΔC₀/2·(1,333)²·e^(-1,333) = ΔC₀·0,5·1,778·0,264 = 0,235·ΔC₀
- Số hạng 2: ΔC₀/2·(1,000)²·e^(-1,000) = ΔC₀·0,5·1,000·0,368 = 0,184·ΔC₀
C3_sốc = (0,235 − 0,184)·ΔC₀ = 0,051·ΔC₀
Với ΔC₀ = 2 × COD_bình thường = 2000 mg/L: C3_tăng thêm = 0,051 × 2000 = 102 mg/L
Hệ thống 3 bể nối tiếp đã “lọc bỏ” phần lớn sốc tải trọng — chỉ 5% ΔC₀ xuất hiện ở đầu ra sau 8 giờ. Đây là khả năng chống sốc tải trọng ưu việt của chuỗi CSTR nối tiếp.
6.2 So Sánh Đáp Ứng Sốc Tải Giữa 1 Bể Và 3 Bể
1 bể CSTR (PT 7-10): Khi C₀ tăng ΔC₀ đột ngột, đầu ra tăng ngay: ΔC1(t) = ΔC₀·(1/τ)/α · (1 − e^(-α·t))
Với α = k + 1/τ và k nhỏ hơn 1/τ: α ≈ 1/τ: ΔC1(t) ≈ ΔC₀·(1 − e^(-t/τ))
Tại t = 2 giờ (H = 2/18 ≈ 0,111): ΔC1 ≈ ΔC₀·(1 − e^(-0,111)) ≈ 0,105·ΔC₀ = 210 mg/L
3 bể CSTR nối tiếp: Tại t = 2 giờ, H = 0,111: ΔC3 ≈ C₀·(nH)²/2! · e^(-nH) = ΔC₀·(0,333)²/2·e^(-0,333) = ΔC₀·0,0555·0,717 = 0,040·ΔC₀ = 80 mg/L
So sánh: Tại t = 2 giờ sau sốc tải: 1 bể → 210 mg/L tăng thêm; 3 bể → 80 mg/L tăng thêm. Chuỗi 3 bể giảm sốc tải trọng gần 2,6 lần so với 1 bể — lợi thế rõ ràng cho các hệ thống XLNT có tải trọng không ổn định.
Kết Luận
Nghiệm PT 7-26 cho bể thứ 2 (C2 = C₀·n·H·e^(-n·H)) và công thức tổng quát PT 7-27 (Ci = C₀/(i-1)!·(n·H)^(i-1)·e^(-n·H)) tiết lộ đặc điểm phân bố nồng độ đặc sắc trong chuỗi CSTR nối tiếp: mỗi bể trung gian có nồng độ đạt cực đại tại H* = (i-1)/n rồi suy giảm, tương ứng với “sóng nồng độ” lan truyền qua chuỗi với độ trễ ngày càng tăng.
Mối liên hệ với phân bố xác suất Erlang không chỉ là sự trùng hợp toán học — nó phản ánh bản chất thống kê của quá trình lưu trú và phản ứng trong chuỗi bể, cho phép dùng lý thuyết xác suất để phân tích hiệu suất xử lý và khả năng chống sốc tải trọng. Kết quả này là nền tảng để so sánh định lượng giữa n bể CSTR nối tiếp, CSTR đơn và PFR — chủ đề của bài viết kế tiếp trong loạt bài này.