Mô phỏng số chân vịt tàu thủy theo phương pháp đa vùng tham chiếu sử dụng openfoam

56 Journal of Transportation Science and Technology, Vol 20, Aug 2016 MÔ PHỎNG SỐ CHÂN VỊT TÀU THỦY THEO PHƯƠNG PHÁP ĐA VÙNG THAM CHIẾU SỬ DỤNG OPENFOAM COMPUTATIONAL APPROACH FOR A MARINE PROPELLER BASED ON MULTI REFERENCE FRAME USING OPENFOAM Phan Quốc Thiện, Bùi Khắc Huy, Lê Tất Hiển, Ngô Khánh Hiếu Bộ môn Kỹ thuật Hàng không, Kỹ thuật Tàu thủy, Trường Đại học Bách khoa Tóm tắt: Hoạt động của chân vịt tàu thủy dựa trên chuyển động quay quanh trục và lực tạo ra nhờ tương tác giữa các lá cánh với dòng nước xung quanh. Bằng phương pháp mô phỏng số sử dụng OpenFOAM, trong bài báo này chúng tôi sẽ nghiên cứu các đặc tính của chân vịt trong quá trình hoạt động ổn định của tàu. Theo đó chuyển động quay của chân vịt sẽ được mô hình hóa bằng phương pháp đa vùng tham chiều (Multi - Referent Frame) trong bài toán tĩnh và mô hình rối hai phương trình k - epsilon sẽ được áp dụng để mô hình hóa chuyển động của lưu chất trong điều kiện rối. Trong bài toán này chúng tôi bỏ qua các yếu tố ảnh hưởng của hiện tượng sủi bọt khi đặt chân vịt trong điều kiện hoạt động ổn định. Các kết quả thu được có độ tin cậy phù hợp với các kết quả dựa trên chuẩn thiết kế của chân vịt. Từ khóa: Mô phỏng số, đặc tính thủy động của chân vịt, lưới cho chân vịt. Abstract: It should be known that interaction between marine propeller and the surrounding water creates the propulsive force. By using computational fluid dynamisc based on OpenFOAM (open source software), this research will focus on the characteristics of a marine propeller. In order to conduct the result, the movement of the propeler is described using Multi - Referent Frame method along with k - epsilon turbulent model in several steady state cases. In these cases, the propeller is considered under various operating condition so that the cavitation is ignored. The obtained results are accurate in comparision of theorical and experimental results. Keywords: Numerical simulation, hydrodynamic properties of ship propeller, mesh generation. 1. Giới thiệu Theo thống kê của Tổng cục Đường Thủy Việt Nam, có khoảng 20% lượng hàng hóa lưu thông nội địa do vận tải đường thủy đảm nhận. Song hành với sự phát triển của giao thông đường thủy là yêu cầu ngày một cao về giảm thiểu chi phí hoạt động nâng cao hiệu suất của tàu thuyền. Điều này dẫn đến việc tối ưu cho quá trình hoạt động của phương tiện là một điều cần thiết. Một trong các yếu tố chính ảnh hưởng đến hoạt động của tàu thuyền là chân vịt. Do đó việc nghiên cứu đặc tính của một chân vịt trong thực tế là một yêu cầu cấp thiết. Tương tự như với các cách thức về mô phỏng cho chong chóng khí của thiết bị bay, việc mô phỏng chân vịt tàu thủy dựa trên bài toán tương tác vật lý giữa bề mặt rắn quay xung quanh một trục và dòng lưu chất bao quanh nó. Bề mặt của chân vịt ở đây được giả định là cứng tuyệt đối. Không có bất cứ biến dạng uốn hoặc nén hay biến dạng cục bộ trên bề mặt. Chân vịt được quay quanh trục của nó sẽ gia tốc dòng nước và đẩy dòng nước về phía sau tạo ra phản lực đưa chân vịt chuyển động về phía trước [1]. Chuyển động của nước là chuyển động của lưu chất thông thường và được mô tả bằng các phương trình bảo toàn trong trong hệ phương trình Navier - Stokes. Trong khi đó chuyển động của chân vịt có thể được mô tả bằng nhiều phương pháp như phương pháp đa vùng tham chiếu (Multi - referent frame), phương pháp đơn vùng tham chiều (Single referent frame) hay phương pháp sử dụng lưới động (Moving dynamic mesh) với các bề mặt lưới trượt lên nhau và với các biên lưới trao đổi dữ liệu với nhau. Trong nghiên cứu này chúng tôi chọn phương pháp đa vùng tham chiếu để mô phỏng đặc tính của chân vịt nhằm đưa ra kết quả nhanh chóng. Có nhiều mô hình rối dạng trung bình Reynolds được sử dụng trong bài toán mô phỏng chân vịt tàu thủy. Do đặc trưng chuyển động phức tạp với độ xoáy lớn nên các mô hình rối sử dụng phương trình như k - epsilon hay k - omega được sử dụng nhiều. Các nghiên cứu của Chang [2], Sanchez - Caja [3] sử dụng mô hình k - epsilon, còn TẠP CHÍ KHOA HỌC CÔNG NGHỆ GIAO THÔNG VẬN TẢI, SỐ 20 - 08/2016 57 Guilmineau [4] sử dụng mô hình rối k - omega SST (một biến thể của mô hình rối k - omega) cho các nghiên cứu liên quan và cho kết quả phù hợp. Bensow [5] sử dụng mô hình “Large Eddy Scale” để nghiên cứu dòng trong trường hợp tương tự nhưng kết quả có độ chính xác thấp. Nghiên này của chúng tôi sử dụng mô hình rối cơ bản k - epsilon kết hợp với mô hình tường (wall model) để đánh giá tương tác tại lớp biên. Hiện tượng sủi bọt trên bề mặt của chân vịt gây ra bởi sự hóa hơi của nước khi áp suất cục bộ thấp hơn áp suất hơi bão hòa. Hiện tượng sủi bọt xuất hiện khi sự chênh áp quá lớn ở bề mặt trên của cánh và xuất hiện nhiều trong quá trình khởi động hay tăng tốc của tàu. Thực tế, mẫu chân vịt nhóm B - Wageningen trong quá trình thực nghiệm đã xét đến việc hạn chế sủi bọt. Do đó, trong nghiên cứu này việc khảo sát tính toán chân vịt nhóm B - Wageningen có thể được xem xét bỏ qua hiện tượng sủi bọt. Mục tiêu chính trong bài báo này là đưa ra một mô hình tính toán số phù hợp với độ tin cậy cho chân vịt tàu thủy nhằm tiết kiệm chi phí với việc đánh giá đặc tính của một chân vịt cụ thể trong điều kiện hoạt động ổn định bằng phần mềm mã nguồn mở OpenFOAM và sử dụng phương pháp chia lưới tự động. Các giá trị về phân bố áp suất, vận tốc và các xoáy nước xung quanh chân vịt sẽ được đánh giá để có một cái nhìn tổng quan về tương tác của chân vịt với môi trường xung quanh. Mặc dù bỏ qua các ảnh hưởng về hiện tượng sủi bọt, chúng tôi cũng sẽ có những phân tích xung quanh vấn đề này trong phần đánh giá kết quả mô phỏng số thu được cho mẫu chân vịt tàu thủy nội địa hiện có ở Bộ môn Tàu thủy, Trường Đại học Bách khoa. 2. Mô hình MRFSimpleFoam Mô hình toán học mô tả bài toán tĩnh mô phỏng chân vịt tàu thủy với phương pháp đa vùng tham chiếu bằng MRFSimpleFoam của OpenFOAM. Trong mô hình này miền mô phỏng được chia làm hai phần tương ứng với hai vùng tham chiếu khác nhau. Phương trình mô tả hiện tượng vật lý gồm hai phương trình: Một là phương trình bảo toàn về khối lượng và một là phương trình động lượng trong hệ phương trình Navier-Stokes: I.U 0  (1)  R I I I p . U U U U         (2) Trong đó các giá trị vận tốc ở các vùng tham chiếu khác nhau ứng với hệ trục tọa độ khác nhau. R IU U r  (3) Khi đó vận tốc quay  sẽ quyết định tính chất của các vùng tham chiếu. Khi mà vận tốc quay này có giá trị bằng không thì bài toán này tương ứng với một bài toán tĩnh thông thường. Vùng quay chứa chân vịt sử dụng hệ tọa độ cục bộ quay quanh trục chân vịt và có cùng vận tốc quay với chân vịt. Vận tốc quay này khác không do đó phương trình bảo toàn động lượng có thêm các thành phần động lượng gây ra bởi gia tốc kéo theo và gia tốc coriolis:  I RU 2 U r     (4) Trong phương trình (2), độ nhớt động học  bao gồm hai thành phần là độ nhớt chảy tầng (l) và độ nhớt rối (t) được suy ra từ mô hình rối hai phương trình k - . 3. Xây dựng mô hình hình học, tạo lưới và thiết lập các điều kiện tính toán Hình 1. Sơ đồ khối với giải thuật SIMPLE trong MRFSimpleFoam với vận tốc tương đối. Hình 1 thể hiện sơ đồ khối với giải thuật SIMPLE trong MRFSimpleFoam đưa ra trong bài viết này với mẫu chân vịt áp để phân tích đặc tính là một chân vịt tàu thủy 58 Journal of Transportation Science and Technology, Vol 20, Aug 2016 nội địa thực tế dựa trên chuẩn thiết kế Wageninen B-series tại Bộ môn Tàu thủy, Trường Đại học Bách khoa [6] (xem hình 2). Hình 2. Số hóa hình học của chân vịt thực tế. (a) Mẫu chân vịt thực tế; (b) Mô hình 3D Các thông số hình học cơ bản của của mẫu chân vịt này: Đường kính chân vịt (D) là 0.4 m; số lá cánh (Z) là 3; tỉ số P/D tại 0.7R là 0.9; tỉ số AE/AO là 0.3. Lưới sẽ được chia cho hình học dựng lại trên máy của chân vịt bằng một module của OpenFOAM gọi là snappyHexMesh [7]. Theo đó, miền lưới bên ngoài là một vùng trụ có chiều dài 3.0 m, đường kính 1.6 m; còn miền lưới bên trong là vùng trụ có chiều dài 0.2 m, đường kính 0.58 m (xem hình 3). snappyHexMesh là một module chia lưới tự động. Để tạo lưới ta cần hiệu chỉnh các thông số giới hạn cho phù hợp về phân bổ mật độ, phân vùng và các tính chất của các đơn vị lưới [7]. Cách làm này nhanh chóng và cho ra các thông số cuối cùng của lưới phù hợp hình học của chân vịt khảo sát. Cụ thể, lưới tạo có 1.76 triệu phần tử lưới, tỉ số max aspect ratio là 6.75, tỉ số max skewness là 3.70. Các thông số lưới được kiểm tra bằng modul checkMesh trong OpenFOAM và thỏa mãn các yêu cầu tính toán. Điều kiện biên cho bài toán mô phỏng được thiết lập theo bảng bên dưới: Điều kiện biên Giá trị Inlet  Velocity-Inlet  Turbulent intensity  Eddy viscosity 0 - 4 (m/s) 10% 10 Outlet  Pressure-Outlet 0 (Pa) Symmetry  Symmetry-Plane Symmetry Chân vịt  Wall  Wall function 62.83 (rad/s) Hình 3. Lưới chia tự động trên OpenFOAM. (a) Mặt cắt toàn miền lưới (b) Mặt cắt miền quay chỉ chứa chân vịt (c) Lưới chia trên bề mặt chân vịt Theo đó, vận tốc dòng vào thay đổi từ 0.0 m/s dến 4.0 m/s để có các giá trị khác nhau của hệ số tiến (J, advanced ratio). Trong khi đó chong chóng quay với vận tốc cố định ở 600 vòng/phút ( 62.83 rad/s). Các biên còn lại được mô tả như miền xa vô cùng. Khoảng cách từ các biên này đến chong chóng đủ xa để giảm tác động đến độ chính xác của kết quả mô phỏng. 4. Phân tích và đánh giá kết quả mô phỏng Lưới được chia làm hai vùng riêng biệt với phần quay và phần đứng yên. Chân vịt được đặt trong vùng quay và áp đặt điều kiện biên tường rắn. Khi quay phân bố vận tốc trên bề mặt chân vịt sẽ càng lớn khi vị trí nó càng xa trục quay và tương tác của bề mặt rắn với dòng nước sẽ đẩy dòng về phía sau. (a) TẠP CHÍ KHOA HỌC CÔNG NGHỆ GIAO THÔNG VẬN TẢI, SỐ 20 - 08/2016 59 (b) Hình 4. Trường dòng qua chân vịt ở J = 0.125. (a) Trên bề mặt hút của chân vịt (b) Trên toàn miền mô phỏng của chân vịt Hình 4 biểu diễn rõ hiện tượng trên chứng minh cho tính đúng đắn của mô hình. Ở đây, chân vịt được mô phỏng trong điều kiện vận tốc dòng nước vào mặt hút là 0.5 m/s ở vòng quay của chân vịt là 600 vòng/phút (tương ứng với hệ số tiến J là 0.125). Lưu ý, trường dòng trong toàn miền lưới như mô tả của hình 4.b là không đối xứng do chân vịt ba lá cánh cũng không đối xứng khi dòng dừng. Có thể nhận thấy từ kết quả mô phỏng thu được, mặc dù phân tích toán học trên hai vùng lưới là khác biệt, tuy nhiên do cả hai đều cùng mô tả một hiện tượng vật lý nên sự tách biệt về lưới không gián đoạn tính liên tục của dòng chảy qua biên phân tách giữa hai vùng lưới. Dựa trên chuẩn thiết kế của chân vịt là nhóm B - Wageningen, các thông số đặc tính hoạt động của chân vịt theo chuẩn có thể được đưa ra làm giá trị tham khảo để đối chiếu với kết quả mô phỏng số của các thông số đặc tính hoạt động của chân vịt khảo sát (xem hình 5) [8]. Cụ thể đối với chân vịt khảo sát: điểm hoạt động có hiệu suất cao nhất trong vùng J từ 0.7 đến 0.75 và đạt giá trị là 0.7 (so với nhóm B - Wageningen là vùng J là 0.8 và đạt giá trị 0.732); kết quả hệ số moment xoắn (KQ) giữa mô phỏng và nhóm B-Wageningen trong vùng J từ 0.125 đến 0.8 hoàn toàn tương đồng với sai số cao nhất chưa đến 0.01%; đối với hệ số lực đẩy (KT) sai số cao nhất ghi nhận được là 8.6% tại J là 0.8. Vùng hoạt động J mà ở đó chân vịt không còn tạo ra lực đẩy ghi nhận được từ mô phỏng là 1.0 so với nhóm B - Wageningen là 1.03. Sự khác biệt về KT và hiệu suất của chân vịt khảo sát so với nhóm B - Wageningen trong vùng J từ 0.8 đến 1.03 có thể được lý giải từ hình học của chân vịt khảo sát không hoàn toàn tương đồng với nhóm B - Wageningen [6]. Liên quan đến hiện tượng sủi bọt, từ các kết quả phân bố áp suất trên mặt hút và trên mặt đẩy của chân vịt khảo sát thu được từ mô phỏng (xem hình 6). Có thể nhận thấy lực đẩy của chân vịt tỉ lệ với độ chênh lệch của áp suất phân bố ở mặt hút và mặt đẩy của các lá cánh. Chênh lệch áp suất càng cao thì lực đẩy sinh ra càng lớn. Mặt đẩy trường dòng bị nén tạo áp suất cao hơn so với môi trường còn mặt hút thì ngược lại tạo áp suất thấp hơn. Và do đó, hiện tượng sủi bọt sẽ xuất hiện nhiều nếu như vùng áp suất thấp ở mặt hút chênh lệch quá lớn vượt ngưỡng hóa hơi. Phân bố trên hình 6.a ghi nhận thấy các vùng có áp suất thấp nhất là các vùng cạnh rìa gần mút của các lá cánh ở mặt hút. Như vậy đây là các vùng sẽ chịu ảnh hưởng trước tiên do hiện tượng sủi bọt. Điều này cũng phù hợp với ghi nhận hiện tượng sủi bọt từ thực nghiệm. Hình 5. Kết quả mô phỏng đặc tính hoạt động của chân vịt khảo sát. 60 Journal of Transportation Science and Technology, Vol 20, Aug 2016 Hình 6. Phân bố áp suất trên bề mặt chân vịt. (a) Áp suất ở mặt hút; (b) Áp suất ở mặt đẩy 5. Kết luận Bài viết đã trình bày việc xây dựng mô hình số phân tích đặc tính hoạt động của chân vịt tàu thủy trên OpenFOAM bằng phương pháp đa vùng tham chiếu với lưới chia tự động. Mô hình mô phỏng số đưa ra khi áp dụng phân tích đặc tính hoạt động của một mẫu chân vịt thực tế đã cho thấy kết quả có độ tin cậy, và có thể hướng đến ứng dụng vào việc phân tích ngược các thông số hoạt động đặc trưng của chân vịt tàu thủy từ hình học của nó. Từ đó, đưa ra các kết quả tham khảo cho việc thiết kế, lựa chọn hệ thống đẩy phù hợp cho phương tiện thủy. Các phân bố vận tốc và áp suất cũng chỉ ra được vùng xuất hiện và một vài ảnh hưởng của hiện tượng sủi bọt lên chân vịt tàu thủy. Bài viết cũng đã giải quyết được các hạn chế về yếu tố thời gian khi đưa ra được một mô hình phân tích phù hợp giúp giải quyết nhanh bài toán mô phỏng số cho chân vịt tàu thủy. Trong các nghiên cứu tiếp theo chúng tôi sẽ tiếp tục phát triển bài toán mô phỏng có xét tới ảnh hưởng của hiện tượng sủi bọt để đưa ra các so sánh tương quan và tăng độ chính xác của bài toán. Lời cảm ơn Công trình được thực hiện tại Trường Đại học Bách khoa, ĐHQG - HCM, thông qua đề tài nghiên cứu cấp ĐHQG loại B năm 2015 (mã số: B2015-20-01)  Tài liệu tham khảo [1] J.S. Carlon, Marine Propeller and Propulsion, Butterworth – Heinemann Ltd., 1994. [2] Chang B., Application of CFD to P4119 propeller, 22nd ITTC Propeller RANS/Panel Method Workshop, France, 1998. [3] Sanchez-Caja, A P4119 RANS calculations at VTT, 22nd ITTC Propeller RANS/Panel Method Workshop, France, 1998. [4] E. Guilmineau, G.B. Deng, A. Leroyer, P. Queutey, M. Visonneau and J. Wackers, Wake simulation of a marine propeller, 11th World Congress on Computational Mechanics, 2014. [5] Rickard E. Bensow and G¨oran Bark, Simulating Cavitating Flows With Les In OpenFOAM, ECCOMAS CFD, 2010. [6] Ngô Khánh Hiếu, Lê Tất Hiển, Đặc trưng hình học và đặc tính thủy động lực chân vịt phương tiện thủy nội địa cỡ nhỏ, Tạp chí Phát triển khoa học và công nghệ, Đại học Quốc gia Tp. HCM, K7-2015, 110-116, 2015. [7] Bùi Khắc Huy, Phan Quốc Thiện, N.K Hieu, Semi-Automatic Meshing Method for Simulating Propeller, Tạp chí Khoa học công nghệ Giao thông vận tải, Đại học Giao thông vận tải Tp. HCM, số 18, 2/2016. [8] Bùi Khắc Huy, Khảo sát đặc tính lực đẩy chân vịt tàu thủy của tàu sông nhỏ, Luận văn thạc sĩ chuyên ngành Kỹ thuật Hàng không, Trường Đại học Bách khoa, ĐHQG Tp. HCM, 01/2016. Ngày nhận bài: 24/06/2016 Ngày chuyển phản biện: 27/06/2016 Ngày hoàn thành sửa bài: 13/07/2016 Ngày chấp nhận đăng: 20/07/2016