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