Sử dụng phương pháp biên nhúng (IBM) xây dựng mô hình 3D tính toán khí động cho cánh máy bay

Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số Kỷ niệm 55 năm Viện KHCNQS, 10 - 2015 101 SỬ DỤNG PHƯƠNG PHÁP BIÊN NHÚNG (IBM) XÂY DỰNG MÔ HÌNH 3D TÍNH TOÁN KHÍ ĐỘNG CHO CÁNH MÁY BAY Mai Xuân Cảnh1*, Lã Hải Dũng2, Nguyễn Trang Minh3, Phạm Minh Vương4 Tóm tắt: Trên cơ sở tính toán khí động cho cánh máy bay với mô hình 2D bằng phương pháp biên nhúng (IBM) đã được công bố [1], nhóm tác giả nghiên cứu xây dựng mô hình tính toán cho cánh với mô hình 3D. Trong bài báo này, phương pháp biên nhúng được tiếp tục cải tiến về mặt phương pháp tính, xây dựng mô hình 3D và thuật giải tính toán song song, để có thể tính được bài toán phức tạp với tốc độ tính nhanh hơn. Chương trình tính đã được thử nghiệm trên nhiều mô hình 3D mà đặc biệt là cánh máy bay mà các tác giả trình bày dưới đây. Kết quả tính toán cho thấy mô hình tính toán xác định được các hệ số khí động, cũng như mô phỏng được cấu trúc dòng chảy trên cánh với độ chính xác cho phép. Từ khóa: Tính toán khí động, Hệ số khí động, Cấu trúc dòng chảy, IBM 1. ĐẶT VẤN ĐỀ Cơ sở lý thuyết của phương pháp biên nhúng và việc áp dụng phương pháp này để nghiên cứu tác động của dòng chảy lên vật thể hai chiều đã được các tác giả công bố [1]. Trong bài nghiên cứu trước, các tác giả đã trình bày các kết quả về vận tốc, đường dòng và phân bố áp suất ở nhiều góc tấn khác nhau. Các kết quả đã khẳng định độ chính xác của phương pháp biên nhúng dùng cho trường hợp dòng chảy qua cánh 2D. Trong phần nghiên cứu này, các tác giả tiếp tục ứng dụng phương pháp tính biên nhúng cho mô hình cánh 3D để có thể tính toán được dòng xoáy, phân bố áp suất cũng như áp suất trên cánh. Phần trình bày cụ thể về phương pháp và các kết quả tính sẽ được giới thiệu dưới đây. 2. GIỚI THIỆU PHƯƠNG PHÁP TÍNH Phương pháp biên nhúng sử dụng cho nghiên cứu này được phát triển từ phần nghiên cứu tính toán khí động với mô hình cánh 2D mà các tác giả đã nghiên cứu trước đó [1]. Điều khác biệt là trong phần nghiên cứu này các tác giả xây dựng mô hình cho bài tính toán 3 chiều. Về mặt lý thuyết, việc áp dụng phương pháp biên nhúng để tính toán khí động cho cánh máy bay nhờ các điểm lực và điểm ảo. Đây là các điểm ở hai phía gần bề mặt. Chia miền tính toán thành hai miền: miền vật lý và miền ảo. Miền vật lý là miền của dòng chảy. Miền ảo nằm bên trong vật thể. Các điểm lực là các điểm của dòng chảy, nhưng không thể có được các giá trị tại các điểm này vì chúng ở gần bề mặt và sự rời rạc là không khả thi. Vì vậy, phương pháp IBM sẽ tính toán các giá trị tại các điểm này mà không rời rạc. Còn các điểm ảo là các điểm trong vật thể và chỉ nhằm mục đích phục vụ tính toán chứ không có nghĩa vật lý. Giá trị tính toán trên các điểm ảo được ngoại suy từ các giá trị trên biên và trên các điểm nằm trên miền lưu chất. Các phương trình tổng quát chú trọng đến sự bảo toàn mô men động lượng và khối lượng. Một thành phần lực trên biên fi được thêm vào phương trình mô men động lượng [4], [6]: 0,j j u x    (1) iji i j Fu f t x      (2) Tên lửa & Thiết bị bay M.X.Cảnh, L.H.Dũng, ,“Sử dụng phương pháp biên nhúng cánh máy bay.” 102 Ở đây: - ui,j (i, j = 1, 2, 3) là thành phần vận tốc theo các trục; - xj (j=1, 2, 3) là tọa độ không gian 3 chiều. - Fij = uiuj + Pδij - ν( j j j i u u x x      ) (3) Trong đó: - P là áp suất dòng chảy có mật độ dòng ρ, độ nhớt động học ν = μ/ ρo; - δij là tensor đường chéo đơn vị (Kronecker delta). Việc giải phương trình mô men động lượng tại mỗi thời điểm phụ thuộc vào điều kiện biên. Lực là một hàm của thời gian phụ thuộc vào vị trí và vận tốc của dòng chảy. Thực tế xi không trùng khớp với các điểm nút, nhưng lực phải được ngoại suy từ các điểm nút này. Lực đặt trên biên nhúng fi bằng 0 trong vùng dòng chảy và khác không ở trong vùng điểm ảo [2], [3], [5]. Lực trên bề mặt để xác định biên vật thể được đặt vào trong miền tính được tính theo biểu thức dưới đây [1], [6]: 1n n i V u f R SH t       (4) Trong đó: - Vn+1 là vận tốc tại điểm lực được nội suy từ vận tốc trên bề mặt của biên; - RSH là thành phần bên phải của phương trình Navier - Stockes trong đó có chứa thành phần đối lưu và lực nhớt; - un là vận tốc tại thời điểm n; Δt là bước thời gian tính toán. Trình tự giải bài toán bằng phương pháp biên nhúng được tiến hành theo các bước sau: Xác định các điểm lực; nội suy để tìm ra các giá trị lực cần thiết để áp đặt điều kiện biên; tính toán vận tốc ở bước dự đoán; giải phương trình Poisson để đảm bảo các điều kiện liên tục; cập nhật trường vận tốc cho các bước tiếp theo. 3. SƠ ĐỒ THUẬT TOÁN Phương pháp giải bài toán theo phương pháp của Mac-Cormack, gồm hai bước tính toán vận tốc dự đoán (n) và tính toán vận tốc chính xác (n+1). Sơ đồ thuật toán tổng quát được thể hiện như trên hình 1. Hình 1. Sơ đồ thuật toán tổng quát. Sau khi xây dựng mô hình tính toán, các điểm tính được phân bố trên các máy tính để tính toán. Hai bước tính toán chính được thực hiện. Bước phỏng đoán dựa trên việc tính toán vận tốc dự đoán từ phương trình động lượng, ước tính lực trên các điểm lực và giải Tính vận tốc dự đoán Ước tính lực trên điểm lực Giải phương trình áp suất (Poisson) Giải ph.trình Poisson (dựa trên v.tốc c.xác) Tính lực trên điểm lực (dựa trên vận tốc chính xác) Xây dựng mô hình 3D Phân bố điểm trên các CPU Hiệu chỉnh vận tốc chính xác K.tra đ/k ổn định Kết thúc Bắt đầu Sai Đúng Bước phỏng đoán Trao đổi dữ liệu song song Bước hiệu chỉnh Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số Kỷ niệm 55 năm Viện KHCNQS, 10 - 2015 103 phương trình áp suất Poisson. Bước hiệu chỉnh được thực hiện nhờ việc trao đổi dữ liệu với bước phỏng đoán để tính toán lại vận tốc chính xác theo sơ đồ Runge - Kutta bậc 2, sau đó tính toán lực và giải phương trình Poisson dựa trên vận tốc chính xác này. Điều kiện ổn định được kiểm tra để kết thúc quá trình tính toán. 4. MÔ HÌNH VÀ KẾT QUẢ TÍNH TOÁN 4.1. Mô hình tính toán Tiến hành xây dựng mô hình số để tính toán khí động cho cánh máy bay có biên dạng kiểu NaCa 0012, với các thông số thiết kế: + Vận tốc, v (m/s): 53 + Độ dài dây cung gốc cánh, bo (m): 2,15 + Độ cao, H (m): 500 + Độ dài dây cung mút cánh, bk (m): 1,08 +Chiều dài nửa sải cánh, l (m): 4,65 + Độ dài dây cung khí động trung bình, ba, (m): 1,64 + Diện tích nửa cánh, S (m2): 7,5 Các điều kiện biên tính toán được trình bày trên hình 1 trong đó vận tốc đầu vào là 53 m/s. Biên trên và biên dưới của miền tính được sử dụng điều kiện biên đối xứng. Ở đầu ra, điều kiện biên áp suất ra được sử dụng. Miền tính được chia lưới đều với số phần tử lưới là: 105106 phần tử; kích thước lưới dx = 0,0204 m, dy = 0,0125 m, dz= 0,027 m (trong mô hình c = 1,62 m). Bài toán được chạy trên hệ thống máy tính song song với 64 CPU và thời gian tính cho một bài toán là 7 ngày để có được kết quả tính toán ổn định. Hình 2. Mô hình tính toán. Hình 3. Hình lưới trên cánh 3D. 4.2. Kết quả tính toán Tiến hành tính toán cho hai trường hợp với 2 góc tấn khác nhau là 2o và 12o. Phương pháp tính toán dòng không dừng đã được sử dụng với hai bước tính khác nhau. Trước tiên dòng khí được nghiên cứu trong một khoảng thời gian nhất định, gọi là thời gian chuyển Tên lửa & Thiết bị bay M.X.Cảnh, L.H.Dũng, ,“Sử dụng phương pháp biên nhúng cánh máy bay.” 104 tiếp để đảm bảo dòng chảy vào được ổn định. Sau đó, bước kế tiếp là chạy để tính giá trị trung bình về lực nâng và lực cản. Phương pháp mô phỏng dòng rối bằng cách tính các cấu trúc lớn dùng mô hình rối LES đã được áp dụng để có thể bắt được các cấu trúc thực của dòng chảy. Phương pháp này cho phép bắt được hầu hết các cấu trúc của dòng và cho kết quả tốt nhất. Tuy nhiên, để có thể bắt được các cấu trúc dòng và tùy từng trường hợp bài toán, chúng tôi đã đạt đến giá trị y+ trên thành nhỏ hơn 10. Với phương pháp mô phỏng rối LES, thuật giải tính toán không dừng là bắt buộc để có thể bắt được các cấu trúc dòng như đã trình bày ở trên. Các kết quả về mô phỏng dòng quanh cánh, vận tốc, áp suất cũng như cường độ xoáy được trình bày dưới đây: Hình 4. Vận tốc ở các góc tấn khác nhau. Hình 5. Phân bố cường độ xoáy. Hình 6. Phân bố áp suất xung quanh biên dạng cánh. Hình 4 thể hiện kết quả phân bố vận tốc trên một mặt cắt 2D qua giữa sải cánh bên phải. Tương tự hình 5 thể hiện phân bố cường độ xoáy và hình 6 thể hiện phân bố áp suất. Kết quả tính cho thấy rằng dòng chảy sau cánh là dòng chảy rối và các cấu trúc dòng sau cánh đã được xé nhỏ thành các xoáy nhỏ khác nhau. Khi góc tới càng lớn thì vùng xoáy Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số Kỷ niệm 55 năm Viện KHCNQS, 10 - 2015 105 sau cánh càng rộng như trong trường hợp góc tới là 12o. Phân bố áp suất trong trường hợp hình 6 cũng thể hiện được rằng áp suất phía dưới cánh và mũi cánh lớn hơn áp suất phía trên cánh. Đây cũng là nguyên lý cơ bản để tạo ra lực nâng và lực cản trên cánh. Hình 7 thể hiện phân bố hệ số áp suất Cp trên toàn bộ bề mặt cánh ở 2 góc tấn khác nhau là 2o và 12o. Ở đây: 21 2 P P C V  Hình 7. Phân bố hệ số áp suất trên cánh. Hình 8 thể hiện phân bố hệ số lực nâng (CL) và lực cản (CD) trên cánh theo thời gian ở hai góc tới khác nhau. Hệ số lực cản trung bình là 0,041 cho trường hợp góc tấn 2° và 0,16 cho trường hợp góc tấn 12°. Hệ số lực nâng trung bình là 0.07 cho trường hợp góc tấn 2° và 0,6 cho trường hợp góc tấn 12°. Các kết quả tính toán được so sánh với một số kết quả nghiên cứu trước đó của của các tác giả trên các hình dạng khác nhau đã khẳng định được tính chính xác của chương trình tính [1], [4] với sai số nhỏ hơn 8,5%. α = 2o α = 12o Hình 8. Đồ thị hệ số lực nâng và lực cản theo thời gian. 5. KẾT LUẬN Nội dung nghiên cứu trong bài báo này đã xây dựng được mô hình tính toán khí động trên cánh máy bay với mô hình 3D bằng phương pháp biên nhúng (IBM). Các bước lập trình cũng như chương trình tính đã được hoàn thiện. Kết quả cho thấy mô hình tính toán cho phép tính được các hệ số khí động, cũng như cấu trúc dòng chảy trên cánh. Các kết quả tính toán sẽ tiếp tục được kiểm chứng bằng mô hình thí nghiệm mà các tác giả sẽ tiến hành trong thời gian tới. Phương pháp tính toán dòng chảy bằng phương pháp biên nhúng, dùng thuật giải tính toán không dừng, trong trường hợp cánh 3D có thể không phải là một trong các phương pháp hiệu quả nhất. Tuy nhiên, phương pháp này có một ưu điểm mà các phương pháp tính truyền thống chưa thực hiện được, đó là với phương pháp biên nhúng, khi cácnh dịch chuyển thì không cần phải chia lại lưới và điều này đảm bảo độ chính xác cũng như hiệu Tên lửa & Thiết bị bay M.X.Cảnh, L.H.Dũng, ,“Sử dụng phương pháp biên nhúng cánh máy bay.” 106 quả tính toán. Phương pháp này sẽ được nghiên cứu phát triển để tính cho cánh có cánh lái và tính phân bố áp suất trên cánh cho các bài toán đàn hồi khí động. TÀI LIỆU THAM KHẢO [1]. Lã Hải Dũng, Nguyễn Trang Minh, Mai Xuân Cảnh, Đào Duy Trường, Phạm Minh Vương (2013), "Xây dựng mô hình số tính toán khí động trên cánh máy bay bằng phương pháp biên nhúng (IBM)", Tuyển tập Hội nghị Khoa học Cơ học thủy khí, Quảng Bình (2013). [2]. Iaccarino, G. and Verzicco, R., "Immersed boundary technique for turbulent flow simulations", Appl. Mech. Rev., Vol. 56, pp. 331–347, 2003. [3]. Minh Vuong Pham, Frédéric Plourde, and Son Doan Kim, "Large-eddy simulation of a pure thermal plume under rotating conditions", Journal of Physics of Fluids, (2007). [4]. M.V. Pham, F. Plourde *, S.K. Doan (2008) "Turbulent heat and mass transfer in sinusoidal wavy channels", International Journal of Heat and Fluid Flow. [5]. Truong D. D, V. M. Pham, N. T. T. Le, "New modeling in the fluid-structure interaction flutter phenomenon analysis", Tuyển tập Hội nghị Cơ học toàn quốc lần thứ IX , Hà Nội, (2012). [6]. Yu-Heng Tseng , Joel H. Ferziger(2003), "A ghost-cell immersed boundary method for flow in complex geometry", Journal of Computational Physics, pp. 593–623. [7]. Nguyễn Chí Công, Phạm Minh Vương, Lê Thị Minh Nghĩa "Mô phỏng dòng qua hình trụ vuông bằng phương pháp LES-Mô hình SMAGORINSKY", Tuyển tập Hội nghị Khoa học Cơ học thủy khí, Hà Tiên (2004). ABSTRACT USING THE IMMERSED BOUNDARY METHOD (IBM) IN ORDER TO BUILD THE 3D CACULATION MODEL FOR THE AERODYNAMIC OF AIRCRAFT WINGS Based on the aerodynamic calculations for wing with 2D models by Immersed Boundary Method (IBM) had been published [1], the authors study to build computational models for the wings with 3D models. In this paper, the immersed boundary method is further improved in terms of methodology, building 3D models and parallel computing algorithms to calculate the complex problems with faster calculated speed. The program has been tested on multiple 3D models which are particularly aircraft wing that the authors present below. The results show that the computational model had identified the aerodynamic coefficients, as well as simulating the flow structure on the wing with allowing precision. Keywords: Aaerodynamic calculations, Aerodynamic coefficients, Flow structure, IBM. Nhận bài ngày 15 tháng 07 năm 2015 Hoàn thiện ngày 05 tháng 9 năm 2015 Chấp nhận đăng ngày 10 tháng 9 năm 2015 Địa chỉ: 1Cục Kỹ thuật Quân chủng PK-KQ; 2 Viện Nghiên cứu phát triển Viettel; 3Viện Khoa học và Công nghệ quân sự; 4Công ty DFM-ENGINEERING. *Email: canhmaixuan75@gmail.com