Mô phỏng số dòng chảy bọt hai pha khí - Nước trong ống đứng

ISBN: 1 Hội nghị Cơ học kỹ thuật toàn quốc Kỷ niệm 35 năm thành lập Viện Cơ học. Hà Nội, 09/04/2014 Mô phỏng số dòng chảy bọt hai pha khí-nước trong ống đứng Dương Ngọc Hải1, 2, Nguyễn Tất Thắng2, Lê Minh Thành2 1 Trường Đại học Công nghệ - Đại học Quốc gia Hà Nội, 144 – Xuân Thủy, Cầu Giấy, Hà Nội 2 Viện Cơ Học – Viện Hàn lâm Khoa học và Công nghệ Việt Nam, 264 – Đội Cấn, Ba Đình, Hà Nội Email:ntthang@imech.ac.vn; thanhlm_vc@yahoo.com Tóm tắt Báo cáo trình bày các kết quả xây dựng mô hình, mô phỏng và tính toán số dòng chảy bọt hai pha khí-nước trong ống đứng sử dụng công cụ phần mềm Ansys FLUENT (version 15). Phân bố của các pha, cấu trúc rối của dòng chảy bọt hai pha khí-nước trong ống có đường kính trong 50mm, chiều dài 3000mm đã được mô phỏng sử dụng mô hình hai chất lỏng (two fluid model) cho dòng chảy nhiều pha. Trong mô hình hai chất lỏng các phương trình bảo toàn khối lượng, động lượng, năng lượng... được giải cho từng pha chất lỏng riêng biệt. Mô hình đóng kín rối k-ε cho từng pha (k-ε dispersed turbulence model) được áp dụng. Profile tỉ phần pha khí (volume fraction or void fraction), vận tốc pha lỏng, pha khí đã được tính toán và so sánh với các kết quả thực nghiệm. Từ khóa: Dòng chảy bọt hai pha, mô hình hai chất lỏng, mô phỏng số, tỷ phần khí, profile vận tốc. 1. Giới thiệu Dòng chảy bọt hai pha có mặt trong nhiều lĩnh vực khoa học, trong các ngành công nghiệp (năng lượng, hóa chất, dầu khí) và trong tự nhiên. Một số ví dụ điển hình bao gồm dòng chảy bọt hai pha quanh các thanh nhiên liệu lò phản ứng nguyên từ nước sôi, trong hệ thống tải nhiệt, sinh hơi của các nhà máy năng lượng, nhà máy điện hạt nhân, trong các cột bọt trong công nghiệp hóa chất, tuyển khoáng v.v. Đặc trưng của dòng chảy bọt hai pha khí/nước (hoặc hơi/lỏng) là sự tồn tại của các bọt khí riêng biệt chuyển động phân tán liên tục cùng pha lỏng. Đường kính lớn nhất của bọt nhỏ hơn nhiều lần so với đường kính ống. Tồn tại các mặt phân cách gián đoạn giữa hai pha. Các mặt phân cách chuyển động và biến dạng liên tục theo không gian và thời gian cùng với các tương tác phức tạp giữa các pha, giữa các bọt và giữa các bọt với chất lỏng [1]. Trên thế giới, nghiên cứu về dòng chảy bọt hai pha đã được tiến hành từ lâu. Cũng như nhiều ngành nghiên cứu khác, các phương pháp mô phỏng số hoặc đo đạc thực nghiệm được áp dụng rộng rãi. Phương pháp mô phỏng số dòng chảy bọt hai pha dựa trên việc giải xấp xỉ các phương trình toán học mô tả dòng bọt. Hiện nay, phương pháp này đã đạt được những tiến bộ đáng kể cùng với sự cải thiện mạnh mẽ khả năng tính toán của máy tính và sự xuất hiện của các phần mềm mô phỏng chuyên dụng dựa trên những kết quả nghiên cứu lý thuyết mới nhất. Tùy thuộc vào sự liên kết giữa các pha, có ba mô hình dòng chảy hai pha gồm: mô hình đồng nhất (hai pha có cùng vận tốc và liên kết mạnh, từ đó hỗn hợp hai pha được giải như dòng một pha), mô hình hỗn hợp (hai pha có vận tốc khác nhau nhưng cũng có liên kết mạnh, hỗn hợp hai pha được giải như dòng một pha có xét đến khác biệt vận tốc giữa hai pha), mô hình hai chất lỏng (mô hình tổng quát nhất, các pha liên kết yếu, hỗn hợp hai pha được giải riêng rẽ và liên kết thông qua tương tác ở biên). Một số công cụ mô phỏng số phổ biến hiện nay đã được áp dụng nghiên cứu dòng bọt hai pha trong ống đứng gồm Ansys CFX (có trang bị mô hình hai chất lỏng, mô hình đồng nhất), Ansys FLUENT (gồm các mô hình hỗn hợp, mô hình hai chất lỏng và mô hình VOF cho dòng nhiều pha dạng phân tách), PHOENICS (chỉ có mô hình hai chất lỏng) [2]. Một số tính toán mô phỏng số dòng chảy bọt hai pha sử dụng các mô hình trên đã được thực hiện [3,4]. Trong nghiên cứu số dòng chảy bọt hai pha, luôn luôn có sự xuất hiện chuyển động tương đối của pha này đối với pha kia do chênh lệch về mật độ và độ nhớt giữa hai pha. Vì vậy hai trường vận tốc riêng biệt cần phải được xem xét trong bất kì nghiên cứu nào của dòng chảy bọt hai pha (mô hình không đồng nhất). Thông thường, dòng chảy bọt hai pha khí-nước có liên kết yếu giữa các pha, mô hình hai chất lỏng cần phải được áp dụng [1,2]. Đối với dòng chảy bọt đẳng nhiệt giữa các pha không có truyền tải nhiệt và năng lượng, mô hình hai chất lỏng có thể được đơn giản hóa với các phương trình bảo toàn khối lượng, động lượng cho từng pha với các đại lượng truyền tải bề mặt giữa các pha. Mục tiêu của nghiên cứu này là thiết lập một mô hình mô phỏng sử dụng phần mềm Ansys FLUENT phục vụ nghiên cứu dòng bọt hai pha khí-nước trong ống đứng. Mục tiêu tiếp theo là kết nối giữa nghiên cứu thực nghiệm và mô phỏng số. Các kết quả mô phỏng số thu được góp phần cung cấp thông tin định tính trong việc thiết kế, phát triển mô hình thí nghiệm, định hướng cho đo đạc. Hơn thế nữa, mô phong số cung cấp các kết quả Dương Ngọc Hải, Nguyễn Tất Thắng và Lê Minh Thành 2 tính toán dự báo cho các điều kiện dòng chảy không thể tiến hành đo đạc thực nghiệm. Các kết quả tính toán mô phỏng của mô hình đem so sánh với các kết quả thực nghiệm và mô phỏng số đã được công bố. Dòng chảy bọt rối hai pha khí-nước trong ống đứng hình trụ có đường kính trong là 50mm, dài 3000mm đã được tính toán và mô phỏng số sử dụng mô hình hai chất lỏng theo quan điểm Euler kết hợp với mô hình đóng kín rối k-ε chuẩn có trong phần mềm Ansys FLUENT. Dòng chảy được mô phỏng với tỷ phần thể tích khí trong dải ~ 2%. Kết quả thu nhận được bao gồm trường vận tốc của các pha lỏng và pha khí, cường độ rối, phân bố áp suất, nhiệt độ, phân bố tỷ phần thể tích của pha khí tại mặt cắt ngang ống bất kỳ. Từ đó ứng xử cuác các tham số của dòng bọt hai pha trong ống đã được nghiên cứu, dự đoán cho các chế độ dòng chảy khác nhau. Trong dải tỷ phần khí được mô phỏng, quá trình liên kết bọt và tách bọt được bỏ qua. Các kết quả thu được đã được so sánh với kết quả đã được công bố của nghiên cứu khác có liên quan [3]. So sánh cho thấy kết quả tính toán sử dụng mô hình mô phỏng số trên nền tảng Ansys FLUENT cho kết quả tương đương với các kết quả thế giới đã công bố. Từ đó mô hình được áp dụng cho các tính toán thử nghiệm và sẽ được đối chiếu với kết quả đo đạc thực nghiệm. 2. Phương pháp số trang bị trong Ansys FLUENT 2.1. Hệ các phương trình Trong tính toán, mô phỏng số sử dụng phần mềm Ansys Fluent, mô hình hai chất lỏng Euler đã được áp dụng. Với mô hình hai chất lỏng Euler các phương trình bảo toàn khối lượng, động lượng được áp dụng giải cho từng pha riêng biệt theo quan điểm Euler. Đối với pha nước được xem như pha liên tục và pha khí như là pha khuếch tán, không có sự trao đổi nhiệt, trao đổi khối lượng trên bề mặt giữa các pha. Hệ các phương trình bao gồm: - Phương trình liên tục đối với pha q:  . 0qq q u    (1) với αq là tỷ phần thể tích, ρq là mật độ khối lượng. - Phương trình cho tỷ phần thể tích của từng pha: 1q p   (2) - Phương trình bảo toàn động lượng đối với pha q:      1 . . n qq q q p qq q q q q pq p u u p g K u u F                      (3) - Phương trình trao đổi động lượng giữa các pha có sử dụng giá trị của hệ số trao đổi trên bề mặt Kpq. Với dòng chảy bọt 2 pha khí-nước Kpq có dạng: 6 p pq p i p f K d A    (4) với dp là đường kính của bọt khí; Ai là mật độ diện tích mặt phân tách (Interfacial Area Concentration - IAC); f là hàm lực cản. Mật độ diện tích mặt phân tách (IAC) là một thông số quan trọng để tính toán quá trình tương tác, trao đổi giữa các pha và được định nghĩa là tỷ số của diện tích bề mặt giữa hai pha trên một đơn vị thể tích của hỗn hợp. Để tính toán mô hình IAC thông thường có hai cách: hoặc sử dụng một biểu thức đại số liên hệ với giá trị đường kính bọt không đổi hoặc sử dụng phương trình truyền tải cho IAC trong đó có xét tới các quá trình liên kết bọt, tách bọt, sự thay đổi kích thước của bọt. Ứng với các vận tốc khí và nước được áp dụng để tính toán trong nghiên cứu này thì có thể được sử dụng một biểu thức đại số liên hệ với một giá trị đường kính bọt không đổi [9]. Đối với pha khuếch tán p, tỷ phần αp, giá trị tham số IAC được tính toán theo mô hình Particle [9]: 6 p i p A d   (5) với αp là tỉ phần của pha khí. Hàm lực cản f được tính toán dựa vào mô hình Schiller-Naumann [9]: 24 DC Ref  (6) Mô phỏng số dòng chảy bọt hai pha khí-nước trong ống đứng 3 với hệ số cản CD và số Re tương đối. Trong đó:  0.68724 1 0.15 / 1000 0.44 1000 D Re Re Re C Re       (7) | |q p q p q u u d Re     (8) với uq, up lần lượt là vận tốc của pha lỏng và pha khí. Đại lượng qF  trong phương trình (3) biểu diễn tổng của các lực trên bề mặt, sự hình thành nên các dạng profile tỉ phần thể tích pha khí phần lớn bị ảnh hưởng bởi các lực này, chúng tác dụng vuông góc với chiều dòng chảy. Đối với dòng chảy bọt trong ống đứng, những bọt khí nhỏ (db<5.8mm) có xu hướng di chuyển hướng về phía thành ống, tuy nhiên những bọt lớn (db>5.8mm) lại có xu hướng di chuyển vào vùng tâm ống. Hiện tượng này được xem xét dưới ảnh hưởng của lực nâng Fl (lift force) được gây nên bởi gradient trường vận tốc chất lỏng theo phương bán kính. Nó tác động vuông góc đối với chiều chuyển động tương đối giữa hai pha.    , ,l q l p q p ql q pF F C u u u             (9) với Cl là hệ số nâng, ρq là mật độ khối lượng pha lỏng, αp là tỷ phần thể tích của pha khí. Đối với vật rắn hình cầu, hệ số Cl = 0.5 đã được tính toán. Đối với dòng chảy bọt hai pha khí nước, thì các bọt khí có hình dạng thay đổi bất thường từ dạng hình cầu, do đó hệ số Cl = 0.05 được lựa chọn [2] Để xét tới vai trò ảnh hưởng của rối trong chuyển động của bọt khí, chúng ta xét lực Ftd (turbulent dispersion force). Lực Ftd được tính toán sử dụng mô hình Lopez de Bertodano (1991) [9] , ,td q td p TD q q pF F C k       (10) với ρq là mật độ của pha lỏng, kq là năng lượng rối động học của pha lỏng, và αp là gradient tỷ phần thể tích của pha khí. Hệ số CTD = 1. Các kết quả thực nghiệm cho thấy các bọt khí tập trung ở vùng gần biên thành ống, nhưng không phải ngay lập tức tiếp giáp với thành ống, gây ra phần tỷ phần thể tích khí thấp ở vùng lân cận thành ống. Theo như Antal (1991) [9] hiện tượng này do ảnh hưởng của lực Fwl (wall lubrication force) nó có xu hướng đẩy các bọt khí ra xa thành ống. Lực này là kết quả của điều kiện không trượt tại biên thành ống. Lưu lượng dòng chảy chất lỏng giữa các bọt khí và thành ống là nhỏ hơn so với dòng chảy giữa bọt khí và phần phía ngoài xa thành ống. Một bọt khí có đường kính db tại một vị trí cách thành ống một khoảng yw chịu tác dụng bởi lực Fwl có dạng: 2 1 2| | , 0 b wl q p wq p w d F u u max C C n y               (11) với |uq - up| là vận tốc tương đối giữa các pha, ρq, αp lần lượt là mật độ của pha lỏng và tỷ phần pha khí. yw là khoảng cách từ bọt khí đến biên gần nhất. Các hằng số C1 = -0.0064, C2 = 0.016 từ kết quả thực nghiệm của Krepper [8]. Lực Fwl = 0 khi yw thỏa mãn điều kiện yw > (Cw2/Cw1)db. 2.2. Mô hình đóng kín rối So với dòng chảy một pha, để mô tả các hiệu ứng rối của dòng hai pha thì số các đại lượng cần được giải trong phương trình động lượng là rất lớn, do đó mô phỏng rối dòng chảy hai pha là vô cùng phức tạp. Vì vậy mô hình rối k-ε dispersed được áp dụng để mô phỏng rối dòng hai pha. Rối của pha lỏng thu được bởi mô hình k-ε chuẩn sử dụng hai phương trình truyền tải năng lượng động học rối k và tốc độ tiêu tán rối ε có bổ sung thêm các đại lượng truyển tải động lượng rối bề mặt. Đối với pha khí các đại lượng rối được tính toán thông qua lý thuyết tương quan Tchen [9]. 3. Thiết lập mô hình mô phỏng số sử dụng phần mềm Ansys Fluent 3.1. Thiết lập mô hình Tính toán được thực hiện ở trạng thái dừng trên một miền tính toán hình trụ ba chiều với đường kính là 50mm, dài 3000mm. Mô hình lưới được tạo bằng công cụ tạo lưới hiện đại Meshing có trong bộ phần mềm Ansys FLUENT. Toàn bộ lưới gồm có 80000 phần tử có dạng hình hộp không bằng nhau được chia đều trong toàn miền (Hình 1a). Độ dày của lớp lưới đầu tiên từ thành ống được chia với một giá trị y+ trong khoảng từ 30 đến 40 để đạt được độ ổn định của bài toán và có được mô hình lực Fwl chính xác (Hình 1b) [9]. Hệ các phương trình được giải quyết trên từng phần tử trong miền tính toán bằng phương pháp thể tích hữu hạn. Phương pháp nội suy upwind được sử dụng để rời rạc hóa các phương trình để thu được các kết quả tốt. Dương Ngọc Hải, Nguyễn Tất Thắng và Lê Minh Thành 4 3.2. Điều kiện biên Đối với điều kiện biên đầu vào của ống, các giá trị vận tốc thực của từng pha cần cung cấp bằng vận tốc superficial của pha đó chia cho tỷ phần thể tích tương ứng, tỷ phần thể tích của pha khí là một hằng số để có được lưu lượng khí mong muốn vào trong ống, đường kính trung bình của bọt d=0.004m; áp suất tại biên ra bằng áp suất khí quyển P=101325Pa (áp suất không khí); biên cứng thành ống là biên tĩnh, không trượt đối với cả hai pha khí và nước. Các đại lượng năng lượng động học rối k và tốc độ tiêu tán rối ε đầu vào được tính toán từ các biểu thức có dạng:   23 2 ink UI (13)   13/4 3/2 0.07in C k D   (14) Với Cμ = 0.09, cường độ rối I=0.16Re-1/8 a, b, Hình 1. Mô hình lưới tính nhìn ngang (a), mặt cắt ngang (b) 4. Một số kết quả tính toán Mô hình đã được kiểm tra thông qua việc tính toán các trường hợp theo như Kepper (1991). Các kết quả tính toán thu được đã được so sánh kiểm nghiệm với kết quả đo đạc thực nghiệm. Trên cơ sở đó, mô hình được sử dụng để tính toán và dự báo tính chất của dòng chảy khi các điều kiện dòng chảy thay đổi, làm tiền đề cho các phương pháp đo đạc thực nghiệm. Hình 2 và Hình 3 biểu diễn phân bố tỷ phần thể tích của pha khí, vận tốc của khí và nước theo phương bán kính của ống ứng với trường hợp Vl=0.4m/s, Vg=0.01m/s, đường kính đồng nhất của bọt khí d=0.004m. a, b, Hình 2. Profile tỷ phần thể tích pha khí tại những vị trí mặt cắt khác nhau của ống a) z = 0.1m, b) z = 3m (Vl = 0.4m/s, Vg = 0.01m/s) Mô phỏng số dòng chảy bọt hai pha khí-nước trong ống đứng 5 Profile tỷ phần thể tích của pha khí đã thu được với dạng đạt đỉnh gần biên “wall-peak” (Hình 2). Các kết quả tính toán được so sánh kiểm chứng với dữ liệu thực nghiệm và tính toán bằng phần mềm CFX của Krepper [2]. Trong Hình 2a profile tỷ phần thể tích pha khí của mô hình tại z = 3m đã được so sánh với dữ liệu thực nghiệm và mô phỏng số CFX. Phân bố tỷ phần khí đạt đỉnh ở vùng gần biên thu được kết quả tương đối chính xác gần kết quả thực nghiệm hơn so với kết quả tính toán sử dụng phần mềm CFX. Khí được đưa vào từ đáy ống, ngay lập tức đỉnh của profile tỷ phần thể tích pha khí ở gần biên được hình thành (Hình 2b với z = 0.1 m) kết quả so sánh với Krepper [2]. Tăng khoảng cách từ đáy ống, profile tỷ phần khí cũng thay đổi, các bọt có xu hướng di chuyển từ vùng tâm ống ra gần thành ống, ở vùng sát biên (R > 0.015) tỷ phần khí tăng lên và ở vùng gần tâm ống thì giảm đi rõ rệt (Hình 2b với z = 3 m). Hình 3 biểu diễn kết quả tính toán profile vận tốc của hai pha khí-nước và contour tỷ phần thể tích của pha khí theo phương bán kính tại vị trí z = 3m. Khác với dòng một pha, profile phân bố vận tốc của pha lỏng dòng hai pha không còn tồn tại phân bố vận tốc lớn nhất ở tâm ống, profile có dạng bằng phẳng ở vùng tâm ống. Kết quả tính toán phân bố vận tốc của mô hình tương đối chính xác với kết quả dữ liệu. a) b) Hình 3. Profile vận tốc của pha nước, khí và contour tỷ phần thể tích của pha khí tại z=3 m (Vl = 0.4m/s, Vg = 0.01m/s) Hình 4 biểu diễn profile tỷ phần thể tích của pha khí cho các chế độ dòng chảy khác nhau ứng với vận tốc pha nước Vl = 1m/s. Hình 4. Các kết quả profile tỷ phần thể tích pha khí ứng với Vl = 1m/s Các kết quả tính toán, dự đoán phân bố tỷ phần thể tích của pha khí ứng với từng trường hợp điều kiện đầu vào khác nhau. Tỷ phần khí lớn nhất ở vùng gần biên đã thu được nhưng ở vùng gần tâm ống thì dạng tỷ phần khí lớn nhất không được hình thành. Lý do có thể tìm thấy ở đây là các hiện tượng tách bọt, kết bọt không được xem Dương Ngọc Hải, Nguyễn Tất Thắng và Lê Minh Thành 6 xét tới trong mô hình. Đối với những chế độ dòng chảy như vậy thì một mô hình mở rộng trong đó có xét tới sự thay đổi kích cỡ bọt do các hiện tượng liên kết bọt, tách bọt phải được quan tâm. 5. Kết luận Mô hình mô phỏng số dòng chảy bọt hai pha khí-nước trong ống đứng sử dụng phần mềm Ansys FLUENT đã được xây dựng. Các tham số quan trọng của dòng chảy bọt hai pha khí nước như tỷ phần pha khí và vận tốc của chất lỏng, chất khí đã được tính toán và so sánh tương đối tốt so với các kết quả mô phỏng số và thực nghiệm đã được công bố. Trên cơ sở đó, mô hình cũng đã được áp dụng để tính toán, dự đoán đối với các chế độ dòng chảy khác nhau. Đối với việc tính toán profile tỷ phần khí dòng chảy bọt hai pha đi lên trong ống đứng, việc xét các lực trên bề mặt giữa hai pha là cần thiết. Mô hình mô phỏng dòng hai pha có trong phần mềm Ansys Fluent có thể mô tả được profile tỷ phần khí với dạng “wall-peak” khi các điều kiện biên của mô hình là phù hợp với các điều kiện của thí nghiệm. Các kết quả nghiên cứu trên chỉ giới hạn trong giả định với một giá trị đường kính bọt không đổi. Một mô hình mở rộng có xét tới sự thay đổi hình dạng, kích thước của bọt do ảnh hưởng của hiện tượng liên kết bọt và tách bọt gây ra sẽ được quan tâm trong các nghiên cứu tiếp theo. Tài liệu tham khảo [1] Ishii, M. and Hibiki, T. (2011). Thermo-Fluid Dynamics of Two-Phase Flow. Springer, second edition. [2] Manninen, M., Taivassalo, V., and Kallio, S. (1996). On the mixture model for multiphase flow. [3] E. Krepper. CFD simulations of a bubbly flow in a vertical pipe, Institute of Safety Research, Annual Report 1999 [4] Agranat, V., Kawaji, M. and Chan, A. M. (2000). Modelling of Gas/Liquid Two-Phase Flows in Vertical Pipes Using PHOENICS. Phoenics Journal of Computational Fluid Dynamics and Its Applications, 13(1), pp. 001-012. [5] Xia Wang, Xiaodong Sun. Three-dimensional simulations of air–water bubbly flows. International Journal of Multiphase Flow 36 (2010) 882–890. [6] Prasser, Horst-Michael. Wire-Mesh Sensors for Two-Phase Flow Investigations. Institute of Safety Research (1999): 23. [7] Drew, D. A. and Lahey, R. T. In Particulate Two-Phase Flow. Butterworth-Heinemann. Boston, MA509–566. 1993 [8] Wellek, R. M., A. K. Agrawal and A. H. P. Skelland, 1996. Shapes of liquid drops moving in liquid media. AIChE J. pp.12:854. [9] Krepper, E., D. Lucas, H. M. Prasser, 2005. On the modelling of bubbly flow in vertical pipes, Nuclear Engineering and Design 235: 597-611. [10] Theory Guide - Ansys Fluent 15.